text:
Random GWAS data generator
I was wondering if there is a tool/script/program that randomly generates GWAS data. The purpose of such a tool would be to use it for educational purposes. So you generate some random .ped, .map and other filetypes. Various distributions like the MAF of all the SNPs, X-Chromosome classificaiton error rates and so on could be set before generating the data. General stuff one can then filter out with plink or similar software to increase the quality of the data.
Then a student would analyze these random files in plink and I already know what the outcome should be, therefore I can see right away if he or she did it right.
Does anyone know if such a program already exists? Or is there a somewhat easy way to program this myself?
Maybe it works if I use existing gene-data from NCBI where I use the given MAF to generate n variations the alleles based on that? Just some random SNPs from dbSNP?

Agreed with the above I would use an existing dataset from NCBIs dbSNP here
Random doesn't work well because HW signal will be meaningless as soon as the Jukes-Cantor distance is exceeded (which it will) and plink needs homologues to generate any signal either in 'reasonable' panmixia or clonality. The result will be panmictic beyond even the most panmictic organism ever known. I don't mind having a simulation of 100% clonalty, because 100% clonality does occur, either through inbreeding or binary division, but this simulation would generate the complete opposite result.


text:
Comparison of phylogeny of SARS2 whole genome vs spike genome
I’ve built two phylogeny at minimum evolution in MEGA-X and extracted these as newick files.
My next step is to compare the two phylogeny to look at if the mutations within SARS2 are more prevalent within the spike, or as a global event across the whole genome. In order to do this I have chosen to use VisualTreeCMP. However, when it comes to choosing how to compare these I am very confused by all the possible options.
Could someone please give me some guidance on the most suitable parameters to give a good image, and comparison between the two phylogeny.
please see the newick format below:
((((((((((((((hCoV-19/Scotland/CVR140/2020|EPI_ISL_425684|2020-03-16:0.00026192,hCoV-19/England/CAMB-73211/2020|EPI_ISL_440287|2020-03-21:0.00026198):0.00000084,hCoV-19/Wuhan/HBCDC-HB-06/2020|EPI_ISL_412982|2020-02-07:0.00000000):0.00000034,hCoV-19/Gambia/NPHL-287/2020|EPI_ISL_561300|2020:0.00000000):0.00000372,hCoV-19/England/NORT-286167/2020|EPI_ISL_573313|2020-07-21:0.00025782):0.00002652,(hCoV-19/South_Africa/KRISP-0543/2020|EPI_ISL_498084|2020-07-12:0.00020599,(hCoV-19/England/ALDP-49EC3C/2020|EPI_ISL_559586|2020-05-20:0.00027741,hCoV-19/Belgium/HCM-0324151/2020|EPI_ISL_420408|2020-03-24:0.00024635):0.00005585):0.00021443):0.00002118,(hCoV-19/Colombia/MAG-INS-87968/2020|EPI_ISL_498152|2020-03-25:0.00015372,((hCoV-19/Australia/VIC8930/2020|EPI_ISL_565286|2020-07-29:0.00000000,hCoV-19/Australia/VIC8929/2020|EPI_ISL_565285|2020-07-28:0.00000000):0.00047159,(hCoV-19/Australia/VIC5819/2020|EPI_ISL_518613|2020-07-20:0.00024688,(hCoV-19/Australia/VIC4715/2020|EPI_ISL_519103|2020:0.00000000,hCoV-19/Australia/VIC2165/2020|EPI_ISL_480772|2020-03-19:0.00000000):0.00001495):0.00004719):0.00036022):0.00009570):0.00004876,((((hCoV-19/Wales/PHWC-16BD74/2020|EPI_ISL_573970|2020-09-15:0.00058675,hCoV-19/Spain/PV-IBV-001568/2020|EPI_ISL_467121|2020-03-13:0.00000000):0.00015901,hCoV-19/England/CAMB-7AAD6/2020|EPI_ISL_433837|2020-04-07:0.00009990):0.00003774,hCoV-19/Greece/HPI-6465/2020|EPI_ISL_451311|2020-03-29:0.00009055):0.00010675,((hCoV-19/USA/CA-SR0097/2020|EPI_ISL_445096|2020-03-13:0.00026193,hCoV-19/Iceland/13/2020|EPI_ISL_417765|2020-02-27:0.00026193):0.00020549,(hCoV-19/England/LIVE-A4CDD/2020|EPI_ISL_500083|2020-03-12:0.00007450,(hCoV-19/Wales/PHWC-16C441/2020|EPI_ISL_574056|2020-09-19:0.00022387,hCoV-19/England/SHEF-10B2F49/2020|EPI_ISL_537266|2020-05-21:0.00003797):0.00018736):0.00008143):0.00010981):0.00001973):0.00006712,(hCoV-19/Italy/LOM-UniMI-L134/2020|EPI_ISL_542130|2020-03-17:0.00064152,hCoV-19/Hong_Kong/HKU-200823-001/2020|EPI_ISL_516798|2020-03-26:0.00015299):0.00016120):0.00001351,(hCoV-19/USA/TX-HMH-1807/2020|EPI_ISL_542855|2020-04-20:0.00010829,hCoV-19/USA/CA-ALSR-1462/2020|EPI_ISL_494618|2020-03-20:0.00015349):0.00017136):0.00008773,(hCoV-19/USA/MA-MGH-00617/2020|EPI_ISL_460393|2020-03-28:0.00022036,(hCoV-19/USA/WA-S28/2020|EPI_ISL_417081|2020-03-02:0.00000726,hCoV-19/Singapore/746/2020|EPI_ISL_500560|2020-07-08:0.00025453):0.00022704):0.00003603):0.00001521,(hCoV-19/Spain/MD-IBV-006433/2020|EPI_ISL_510145|2020-04-03:0.00030013,hCoV-19/Austria/CeMM0009/2020|EPI_ISL_419662|2020-03-14:0.00000000):0.00025838):0.00001169,hCoV-19/Wuhan/IME-WH01/2020|EPI_ISL_529213|2019-12-30:0.00000000):0.00000016,hCoV-19/Wuhan/IPBCAMS-WH-01/2019|EPI_ISL_402123|2019-12-24:0.00000000):0.00000045,hCoV-19/USA/WA-S106/2020|EPI_ISL_417159|2020-02-29:0.00026265,NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome:0.00000000);

(((((((((((hCoV-19/England/ALDP-49EC3C/2020|EPI_ISL_559586|2020-05-20:0.00013635,hCoV-19/Belgium/HCM-0324151/2020|EPI_ISL_420408|2020-03-24:0.00006529):0.00000773,hCoV-19/South_Africa/KRISP-0543/2020|EPI_ISL_498084|2020-07-12:0.00026099):0.00013025,(hCoV-19/England/LIVE-A4CDD/2020|EPI_ISL_500083|2020-03-12:0.00016259,(hCoV-19/Gambia/NPHL-287/2020|EPI_ISL_561300|2020:0.00016807,hCoV-19/England/CAMB-73211/2020|EPI_ISL_440287|2020-03-21:0.00013706):0.00000534):0.00000008):0.00000278,(hCoV-19/England/NORT-286167/2020|EPI_ISL_573313|2020-07-21:0.00036517,((hCoV-19/Australia/VIC8930/2020|EPI_ISL_565286|2020-07-29:0.00000000,hCoV-19/Australia/VIC8929/2020|EPI_ISL_565285|2020-07-28:0.00000012):0.00009405,(hCoV-19/Australia/VIC5819/2020|EPI_ISL_518613|2020-07-20:0.00006516,(hCoV-19/Australia/VIC4715/2020|EPI_ISL_519103|2020:0.00003344,hCoV-19/Australia/VIC2165/2020|EPI_ISL_480772|2020-03-19:0.00003367):0.00000193):0.00000601):0.00017105):0.00002318):0.00000716,hCoV-19/Greece/HPI-6465/2020|EPI_ISL_451311|2020-03-29:0.00004599):0.00008629,(hCoV-19/England/SHEF-10B2F49/2020|EPI_ISL_537266|2020-05-21:0.00015654,(hCoV-19/Wales/PHWC-16C441/2020|EPI_ISL_574056|2020-09-19:0.00019758,hCoV-19/Wales/PHWC-16BD74/2020|EPI_ISL_573970|2020-09-15:0.00023814):0.00021545):0.00003087):0.00000571,((hCoV-19/USA/TX-HMH-1807/2020|EPI_ISL_542855|2020-04-20:0.00023928,hCoV-19/Colombia/MAG-INS-87968/2020|EPI_ISL_498152|2020-03-25:0.00013017):0.00003266,(hCoV-19/Scotland/CVR140/2020|EPI_ISL_425684|2020-03-16:0.00003644,hCoV-19/Italy/LOM-UniMI-L134/2020|EPI_ISL_542130|2020-03-17:0.00023391):0.00000153):0.00000121):0.00012657,(hCoV-19/Singapore/746/2020|EPI_ISL_500560|2020-07-08:0.00034258,(hCoV-19/Iceland/13/2020|EPI_ISL_417765|2020-02-27:0.00012787,(hCoV-19/England/CAMB-7AAD6/2020|EPI_ISL_433837|2020-04-07:0.00012752,(hCoV-19/USA/CA-ALSR-1462/2020|EPI_ISL_494618|2020-03-20:0.00010172,hCoV-19/Hong_Kong/HKU-200823-001/2020|EPI_ISL_516798|2020-03-26:0.00020057):0.00000680):0.00004213):0.00002475):0.00002767):0.00000390,(hCoV-19/Wuhan/HBCDC-HB-06/2020|EPI_ISL_412982|2020-02-07:0.00007999,(hCoV-19/Spain/PV-IBV-001568/2020|EPI_ISL_467121|2020-03-13:0.00019054,((hCoV-19/Wuhan/IME-WH01/2020|EPI_ISL_529213|2019-12-30:0.00003592,hCoV-19/Spain/MD-IBV-006433/2020|EPI_ISL_510145|2020-04-03:0.00020087):0.00000317,(hCoV-19/USA/CA-SR0097/2020|EPI_ISL_445096|2020-03-13:0.00011733,(hCoV-19/USA/MA-MGH-00617/2020|EPI_ISL_460393|2020-03-28:0.00013271,(hCoV-19/USA/WA-S28/2020|EPI_ISL_417081|2020-03-02:0.00003240,hCoV-19/USA/WA-S106/2020|EPI_ISL_417159|2020-02-29:0.00006806):0.00001270):0.00000336):0.00008627):0.00001038):0.00001150):0.00004587):0.00000335,hCoV-19/Austria/CeMM0009/2020|EPI_ISL_419662|2020-03-14:0.00013033):0.00000452,hCoV-19/Wuhan/IPBCAMS-WH-01/2019|EPI_ISL_402123|2019-12-24:0.00010072,NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome:0.00000000);

And a link to the tree viewer:
https://eti.pg.edu.pl/treecmp/WEB

I am just requesting the basic parameters, on how one would get the most out of the programme... or failing that another programme with less in depth parameterization.

These are also rooted on the Wuhan reference genome

Interesting, will be at least until the weekend before I can respond

You need a variety of metrics in context, but like ... hmmmm, you're expecting a lot from a single tree comparison tool, albeit its parameter rich.

This is for an undergraduate dissertation, so I mean it doesn’t need to be crazy in depth analysis. Due to the pandemic, we’ve been forced away from wet lab work. So my main issue is having to self teach so much new analytic methods. So just the basic parameters, on how one would get the most out of the programme... or failing that another programme with less in depth parameters

Okay you've just made my "when I get around to doing it" that bit more "when I get around to doing it". Still its an interesting approach and it is a topical paper for sure.

I don’t know whether I should take that as positive or negative? :D

Hi @TheoJones I'll do this over next weekend, apologies for the delay.

source: Stack Exchange

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

sort: votes


text:
How do you calculate the top L/5 score?
I am currently investigating how to predict protein structure and contact map predictions. There, I see things like top L/5 and top L/2 scores as a way to evaluate the contact map. What formulas do you use to evaluate these? Also, what is a meaningful evaluation method?

L is protein length so you divide that by two or five!
Sorry, I could not resist that.
Anyway, there are a few metrics to rank the contact, so depends on the paper. The original EVFold paper is a good place to start https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0028766 even if you are interested in the AlphaFold stuff.

@MatteoFerla Thanks for answering my question, I understand that L/5 is the length of the protein sequence divided by 5. Does that mean that Top L/5 represents the best prediction for all patches of L/5?

@MatteoFerla Also, what does Short, Medium and Long mean in the accuracy column of the forecast? This is the paper: https://academic.oup.com/bioinformatics/article/35/22/4647/5487385

That sounds like the range cutoffs —A range contact is not the distance of atoms, but the distance in residues in the linear sequence —in that paper <12 is short, >23 long. Honestly, I've not read that paper —the Zhang group have made I-Tasser so are a good group. But this is paper dependent, in the AlphaFold one they divide the range in 64 bins, but still refer to short-range/medium-range/long-range.

A contact map predictor will usually output the probability of a contact for every residue pair in a protein, often ommitting very close pairs (e.g. less than 6 residue separation).
To calculate the L/n score, rank the predicted contacts descending by probability and look at the first L/n contacts. The fraction of these that are actually contacts in the native structure is the L/n contact precision score.
Often, you would filter the list to look at only contacts of a certain sequence separation, e.g. long range contacts between residues of >23 separation (definitions vary between papers). If you have a good predictor this score is high for L/10 and falls as you go to L/1, which is expected as proteins don't usually have that many long range contacts to find.

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


text:
How to concatenate two fasta sequences by py
I have three backbone and i want to concatenate 70 sequence into these backbone.
such like:
fasta file 1:
Backbone1: ACTGCGTATGTCAGCATATAGACCGCCT
Backbone2: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTC
Backbone3: GTCATGCTACGTGCTACGTGCTCGTCGCGC

fasta file 2:
sequence_1:CGCTAGTGTAGCTGCTACTCCG

output sequence:
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCGs

Use Biopython's seqio

Can you give me some tips for writing codes?

First of all, those are not fasta sequences. Fasta is a specific format, and requires a > followed by a sequence identifier on one line and the one or more lines of sequence.
Assuming you really have two files in the format you show, you can do:
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]; next} for(b in bbone){print b$2}}' file1 file2
GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

The gsub(/ /,""); just removes the extra space you have in file1 after each :. If you also want to have the sequence names as you show, use:
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]="Backbone"NR; next} for(b in bbone){print bbone[b]"+"$1":", b$2}}' file1 file2
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

If you really need to do this in python for some reason, please edit your question and show us your python script so we have something to work on and expand. I really don't want to write a whole python program for something so trivial.


text:
If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?
I am reading a book about RNAseq analysis and it says

"To calculate the probability that a read will map to a specific gene, we can assume an average gene size of 4000 nt (100 M nt divided by 25,000 genes). At 30 M reads equivalent to 30× coverage, at single read 100 nt (or paired-end read 50 nt) length, we can expect a single read to map to the average expressed and length gene 4000 nt× 30 coverage/100 nt 1200 times. Thus, if the gene is expressed at a level of 1/1200 compared to the average gene, then we have a 50:50 probability to have a read map to it."

I understand that a single read would map to the average expressed 1200 times but I don't understand how the 50:50 probability is calculated. Please could you explain how this probability is determined.

I actually disagree with that.. I guess I write this down as a discussion.
The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.
However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:
plot(dpois(0:10,lambda=1),ylab="Probability",xlab="No of reads")

The probabilities are like:
data.frame(n=0:5,p=dpois(0:5,lambda=1))

  n           p
1 0 0.367879441
2 1 0.367879441
3 2 0.183939721
4 3 0.061313240
5 4 0.015328310
6 5 0.003065662

The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206

I guess the author just assumes uniform distribution instead of poisson.

possible.. but why 50:50 ?

I was curious about this question but didn't know what the distribution would be, so I left a very clumsy monte-carlo simulation running yesterday.  After about 8,500 iterations I get exactly these values.

One other thing: Could the author have taken the probability determined for the lambda=1200 value of getting >=1200 (about 50%) and incorrectly assumed it'd stay the same for lambda=1?

There are a lot of back-of-the-envelope calculations in the original statement. I don't think 50:50 is supposed to be an exact number. The point is that you cannot be very confident in detecting that gene.


text:
Survival analysis using CoxPH - Effect of covariates
Hi and sorry for the long post in advance,
I'm doing a survival analysis of lung cancer patients using Python's lifelines package.
According to the documentation, the function plot_partial_effects_on_outcome() plots the effect of a covariate on the observer's survival. I do understand that the CoxPH model assumes that the log-hazard of an individual is modelled by a  linear function of their covariates, however, in some cases the effect of these covariates seems to perfect when plotted. I'll illustrate these with three images below describing the issues I have in my data analysis.
Example 1:
This one I'm providing as it looks realistic to me, there is a clear difference in the effect of the covariate when positive or negative.

Example 2:
This one just looks gimmicky, my concern is firstly there seems to be a perfect relation between the pathologic stage of cancer and the survival, it's logical that higher stages should impact the survival negatively, however, I don't understand why the graph is so symmetric as the survivability of cancer shouldn't be proportional to the patient cancer staging.

Example 3:
Same problems as for example 2, this time with a binary not categorical variable (so I'm less reluctant on disbelieving it), however, the curves do seem to symmetric, almost like it's the same curve but with a standard deviation. I also highlighted a region of interest which I don't understand why both curves stabilize (maybe lack of patients in those regions?)

I would post to cross-validation unless @StupidWolf has done survival analysis. I've never done it, so can't help. Those are very different probabilities. and the drop is simply the tail of the distribution and sampling at event 5000. I wouldn't worry about that- the difference between the top and bottom graph is ... waoooo

Thank you for your contribution! I'll wait a bit before moving it to cross-validation, just thought that considering the domain it would be more fit here.

You were right to do so. I retract my original note.

How come "condition 1" (red line) in the first plot does not start from 1.0? Can it just be the result of this curve being to steep?

That's actually a good observation, and I'm not sure why it's happening. I've looked into it and zoomed the graph, and it's actually not starting at 1, but at 0.9.
It may be something to do with the variable not being fit for this type of analysis as there is only one positive case in 496. The only positive case is a patient who died on day 24, so I can see the curve being steep however I have no idea why it doesn't start at 1

The actual hazard ratio of the covariate ´paper_Selected.Mutation.Summary.TP53 p.R280I´ is 147.00275 which is actually the highest for the LUSC analysis. Now my only assumption why this might be correct while it seems strange is that there are actually some patients dying after 0 days, although they do not have this covariate positive.
Then again, it might just mean that with only 1 positive case in 496 samples this covariate can act more like an outlier than anything. If anyone has a supposition why it doesn't start at 1 I can provide more info.

Example 2 and 3 seem symmetric because they are forced to by the model. The Cox proportional hazard relies on the proportional hazards assumption. This basically means that the ratio of the hazards are constant over time. If you were to plot the hazard ratio, it would look like a flat line over all Event_at. Another way of saying this is that the hazard ratios are constant over time.
How that shows up in a survival plot is a little different because of the relationship between hazards and survival. The survival curves end up following the pattern exemplified in Example 2.
As for example 3, I would guess that there are few individuals with event times past 3500. This is why you have those large jumps.
Proportional Hazards
So what does the proportional hazards assumption mean? Below are some plots of hazard functions to help your intuition on what the proportional hazards assumption means. In all of the following scenarios, the hazard ratio is 2 for blue-to-red. However, the proportional hazards assumption is only true in A, B, and C.

Since the proportion is remaining constant, it can look a little hard to see. To make it easier on the eyes, you can log-transform the hazard functions. Below is a plot that does that.

As you can see, the hazards are a constant distance from each other for A, B, and C. Basically, the proportional hazards assumption is saying that your scenario looks something like one of those plots.
In scenario D, the hazard ratio varies over time. At the beginning it is basically HR=1. At the end it is basically HR=4. With the proportional hazards assumption, Cox regression will essentially take the 'average' of the hazards over time. In this case it would give you an estimate of HR=2. Depending on the scenario, this 'smoothing' may be more or less acceptable. But this is basically what the proportional hazards assumption does in the background.
The proportional hazards assumption is actually more problematic in the cases where the two hazard functions cross each other. Below is an example of that.

In this scenario, blue has an increased hazard at the start (relative to red) then a decreased at the end. A Cox proportional hazards model would give you HR=1 in this case, because the average over all time points is no relationship.
So how do you tell if you violate the proportional hazards assumption? You may want to use a stratified Kaplan-Meier by your variables to get an idea of what the hazards look like without stipulating the proportional hazards assumption. You can use the Kaplan-Meier estimator to get a non-parametric estimate of the discrete hazard (i.e. it does not make the proportional hazards assumption). Then you can log-transform those hazards to get a plot like the second. There are also statistical tests and other diagnostics (Schoenfeld residuals) to check this assumption.
The next question is do those violations matter? This somewhat depends on the question you are trying to answer and how substantial the violations are. In cases where the two log-hazard curves don't cross (like in D), it may be acceptable to smooth over the variations. However, it settings where the curves cross in the final figure, it is probably best to not ignore this problem.
The original write-up of this answer was here

Please upvote this @pzvich this is a truely excellent response and helps me understand survival analysis


text:
Does BLASTClust guarantee that proteins in different clusters are dissimilar?
I need to find dissimilar proteins. Looking through the PDB I found the weekly BLASTClust results of proteins that are 30% similar. However, I do not know if protein A in cluster 1 is guaranteed to be less than 30% similar to protein B in cluster 2.
The best and only documentation I can find for BLASTClust is here (link).
It states:

The program begins with pairwise matches and places a sequence in a cluster if the sequence matches at least one sequence already in the cluster.

That's good, but consider a set of three proteins (A, B, and C). A is 40% similar to B and 40% similar to C. B is 20% similar to C. Say that the algorithm places B in cluster 1 and then C in cluster 2 as they are dissimilar. What happens when it tries to place A? Because A is similar to both B and C will they be joined into one cluster (this would be the best result for me) or is A just placed into one cluster even though it is similar to a protein in another cluster?

what do you mean "dissimilar"? Do you mean "unrelated and non-homologous" or do you mean "homologous but highly diverged"?

see also here for more documentation: https://extras.csc.fi/blast/doc/blastclust.html

@MaximilianPress Sequence identity is the measure I'm using.

sure, identify is a measure you can use, but it doesn't tell us what question you want to ask with these proteins. pick any two random proteins from uniprot, and they will be dissimilar 99.99+% of the time, because they are simply unrelated. why not do that?

Based on the documentation (that I linked to in comment), it seems that the clustering algorithm used by blastclust is single-linkage hierarchical clustering.
SO what you can expect is that A will be placed with the protein to which it is most similar. And then the last protein will be added as an outgroup to that cluster due to its similarity to A (assuming that there are no other proteins).
So if your goal is "I want clusters that are guaranteed to contain dissimilar proteins", then you might not be able to naively take blastclust output and interpret it as such. You could either postprocess it to find the set of clusters in which no protein is more than 30% similar to a protein in any other cluster (which may just be one big cluster!), or find another clustering algorithm that is more suited to your task.
Update
It does look like blastclust uses a minimum similarity (not identity) threshold (from blastclust docs). So they use an implementation of single-linkage clustering that stops before merging everything into a single cluster:
If the coverage is above a certain threshold
 AND
the score density is above a certain threshold,

these two sequences are considered to be neighbored.

Thus determined neighbor relationships is considered symmetric and provides
the base for clustering by a single-linkage method (which puts a sequence
to a cluster if the sequence is a neighbor to at least one sequence in the
cluster)."

I don't think that this is a general property of the generic clustering algorithm though (see the "Naive algorithm" section).

From your link: "In single-linkage clustering, the distance between two clusters is determined by a single pair of elements: those two elements (one in each cluster) that are closest to each other. The shortest of these pairwise distances that remain at any step causes the two clusters whose elements are involved to be merged." So it looks like all three proteins in my example will be in the same cluster at the end. This means that any proteins in another cluster will be guaranteed to be less than 30% similar to my example cluster.

@lazer-guided-lazerbeam i'm not sure that i follow that argument. yes, all of your proteins will eventually be joined into one "cluster" by hierarchical clustering- the output of a classic hclust algorithm is not a set of clusters but a tree, where you can "cut" at different heights to select different intensities of clustering. so you can just select the root, and it's all one cluster! blastclust presumably implements some heuristic to choose where to cut (silhouette analysis is one such approach). i am not sure that you can guarantee the 30% threshold outside of artificial examples.

The single linkage article states that two clusters will only be merged if their two closest elements are are close enough (in this case more than 30% similar). Therefore every element in one of two separate clusters must be less than 30% similar to every element in the other cluster, right?

@lazer-guided-lazerbeam I don't see that in the wiki article, quite possible I am missing something. I am specifically looking at the "Naive Algorithm" section and the worked example. I do actually see that something similar is in the blastclust docs- which argues that you are right.


text:
Measure the purifying selection for certain taxa along a phylogeny
I am posting this message because I need clarity about the analysis I want to perform.
In my analysis, I have a homologous gene in 13 species, and I would like to evaluate the selection pressure of this gene (its dN/dS) in 3 species of the tree and test if this value is significantly different from the neutral model (dN/dS = 1).
So here is the protocol I am considering, and I would like you to tell me if I am right or not:
To calculate this selection pressure, the ratio of non-synonymous to synonymous substitution rates (dN/dS) must be calculated.
For this I use CODEML implemented in ete3.
So I need :
A phylogenetic tree of the gene (performed with the AA alignment of the gene via Iqtree)

A codon alignment (performed with the DNA alignment of the gene with pal2nal)

Mark branches in the tree.

Here is the exemple where I want to get the dN/dS value of SP5,6 and 7.:
         /-SP1
        /-|
       |   \-SP2
       |
     /-|      /-SP3
    |  |   /-|
    |  |  |   \-SP4
    |   \-|
    |     |   /-SP5
  /-|      \- |  /-SP6
 |  |         \-|
 |  |           |
 |  |            \-SP7
 |  |
 |  |   /-SP8
-|   \-|
 |      \-SP9
 |
 |     /-SP10
 |  /-|
 | |   \-SP11
 |-|
   |   /-SP12
    \-|
       \-SP13

And I want to calculate the dN/dS of SP5,SP6 and SP7 in order to confirme the hypotesis that the gene is under strong purifying selection on this 3 species.
I mark the SP5, SP6 and SP7, run a free-branch model and look at the Codeml output (the dN/dS value should be the same for the 3 marked species right?) Then this value is kind of the mean dN/dS value between three leafs?
Suppose I get an omega = 0.15
If I want to test that this value is significantly different from 1, I need to mark the same sheets and run the neutral model, then compare the two models (neutral vs. free) using chi2, and if the p-value is <0.5, then the dN/dS value (0.15) is significantly different from 1.
I'm right  ?

Hi @chippycentra ... thats a funny username mixing.crossing cultural colloquialisms from very different regions, Ireland and the Northern GB! Anyway your question is fine but you need to provide some code. I get what you are doing, this is Python. Otherwise the response will be equally as vague. I do not think codeml is implemented in ETE3 ... I may have missed recent developments, but I'd be surprised if ETE3's developers have gone down that path.

Hello, ahah I did not know that it was a random username ^.^. Well I used this protocol http://etetoolkit.org/docs/latest/tutorial/tutorial_adaptation.html#test-on-branches) ete3 implemented codeml (they still use codeml but with a python code).  I edited the post in order to be more precise.

There is no way that is a random username ... anyway you are right about CODEML. I want to read how they have done this first because it would imply CODEML (UCL, UK) have teamed up with the ETE3 guys (Spain).  Its possible CODEML code was just given to ETE3 and has changed its user license. I can give you the answer in terms of CODEML outside ETE3, but inside it ... this is new to me.

Well, since ete3 simply write a script for codeml and then run this script on it, you can give my an answer with codeml and i will understand :-)

Okay, not till the weekend, its just that Biopython was the standard wrapper for this sort of stuff so I'm astonished by the ETE3 thing.

it is juste that ete3 allows to deal easily with trees I guess

In my opinion .. it is very good for handling and traversing trees and great in data pipelines. As a CODEML wrapper I never knew it existed.

@Michael Just as an aside, do you know a software that allows, from an alignment of n sequences, to calculate the dN, dS, dN/dS between each pairwise sequence? But with a robust method such as likelihood (and not a counting method).

Please see this post https://bioinformatics.stackexchange.com/questions/14801/multiple-alignment-cost-application/14813#14813 ... I'm not a fan of what I think you are proposing, but I might be wrong in what you are asking

The following is not the clearest explanation of dN/dS humankind has seen, but I hope you get the gist and core idea of likelihood is absolutely correct. I am wondering if you are thinking of a different test to the way I would do this calculation. This ain't a text book clarity answer.
Anyway what you are proposing is a likelihood ratio test (LRT) for  selection based on the parameterisation of omega. Your idea isn't quite correct ... its not dN/dS value you are looking at to assess significance, its the likelihood. This point is critical and is the most important difference in our thinking on this test.
LRT works because its a log likelihood and subtracting two log likelihood is the same as dividing them (I forget the theory - but you subtract the likelihood to work out the p-value via Chi squared). There's some theory about why the Chi-squared distribution works but I can't remember it.
You need to work out the degrees of freedom which is the parameter differences between 'free branch model' and the neutral model and this value is LOW. If the neutral model is the "global model", then its 2 degrees of freedom - IN YOUR TREE. This is because the SP you are looking at form a single clade (monophyly). From memory you stick a #1 on the branch leading into the clade, ie. defining the node and thats your 'model' on the treefile. Here you have two dN/dS one for the whole tree excluding your clade and one for your taxa. This gives you a likelihood.
Note you could stick a #1 on each branch within the clade, i.e. the treefile. That means something different - like every branch was subject to ... in this case purifying selection ... not just the common ancestor.
You now work out the likelihood of the 'global' model... i.e. remove the #1 and make sure it gives you only 1 dN/dS value.
Subtract the likelihoods and multiple by two - if its LRT (super important).
You now need to work out the degrees of freedom. If its the same model but one has an additional node with a different dN/dS ... the df = 2. If you've dumped a #1 on every branch in the clade ... I still think df = 2 - you need to check.
If the difference in likelihood results in p < 0.05 you are correct the difference is significant ... whatever that difference might be.
If the neutral model is some different calculation of omega within the population genetics process (see here) then it will be the parameter differences between your model and the 'neutral' model.
To reiterate
I like simulations so much I type in reiterations ... So again the LRT usual will usually subtract one likelihood from another and THEM multiply the probability by two ... and compare the result on the Chi squared table for the difference in degrees of freedom. And again what you are looking at is not the value of dN/dS but the difference in LIKELIHOOD. That is very important. You need to be careful with the degrees of freedom on an LRT because they can introduce an additional subtraction constant, test dependings (its ages since I've done one so I can't remember).
If the test is significant and one dN/dS is x for the clade and y for the rest of the tree, then that differences is significant and if y is >1 (bit unlikely) and x < 1 then sure thats purifying selection. Usually clades evolve via positive selection rather than purifying selection.
From memory that is my understanding ... but if you have a neutral model of omega and a purifying model of omega ... yeah sort of, but usually thats for a 'global' model, i.e. no differences within the tree. In your case your interested in a specific clade. The idea is still the same its measured through likelihood.
Its so long since I've done this, but I do caution not to confuse two different tests.

I do understand Yang et al did an MCMC for dN/dS ... what I'm thinking of here is maximum likelihood. You'd need to watch for that carefully and thats what you are referring to when you say Chi-squared, because that is LRT 100%. LRT is maximum likelihood.

Thank you very much for your time and all the explanation it wad really usefull !!! I also enjoy it since you know PAML implemented in biopython, I'm currently doing a dN/dS pairwise analysis (runmode=-2), but it doesn't seem to work, do you have an idea ? I've created a post here with more informations  : https://www.biostars.org/p/474752/

Hi @chippycentra much appreciated. I'm afraid I've never done Biostars - I don't even have an account - ... but I will have a look and post my thoughts here. PAML does work but can be a bit fiddly.


text:
Chip-seq down-sampling problem
Control sample chip read stat
56234445 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
55173337 + 0 mapped (98.11% : N/A)
56234445 + 0 paired in sequencing
28078024 + 0 read1
28156421 + 0 read2
53808040 + 0 properly paired (95.69% : N/A)
54626508 + 0 with itself and mate mapped
546829 + 0 singletons (0.97% : N/A)
162718 + 0 with mate mapped to a different chr
112602 + 0 with mate mapped to a different chr (mapQ>=5)

IgG sample chip seq read stat
32745255 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31546173 + 0 mapped (96.34% : N/A)
32745255 + 0 paired in sequencing
16362391 + 0 read1
16382864 + 0 read2
30934194 + 0 properly paired (94.47% : N/A)
31269484 + 0 with itself and mate mapped
276689 + 0 singletons (0.84% : N/A)
79786 + 0 with mate mapped to a different chr
49658 + 0 with mate mapped to a different chr (mapQ>=5)

knockout sample Rep 1 stat
28044731 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
27459222 + 0 mapped (97.91% : N/A)
28044731 + 0 paired in sequencing
14014483 + 0 read1
14030248 + 0 read2
26685254 + 0 properly paired (95.15% : N/A)
27162846 + 0 with itself and mate mapped
296376 + 0 singletons (1.06% : N/A)
116684 + 0 with mate mapped to a different chr
88965 + 0 with mate mapped to a different chr (mapQ>=5)

Knockout sample rep 2 stat
32410094 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
31774629 + 0 mapped (98.04% : N/A)
32410094 + 0 paired in sequencing
16196814 + 0 read1
16213280 + 0 read2
30802808 + 0 properly paired (95.04% : N/A)
31561218 + 0 with itself and mate mapped
213411 + 0 singletons (0.66% : N/A)
89688 + 0 with mate mapped to a different chr
63838 + 0 with mate mapped to a different chr (mapQ>=5)

Do i need to down-sample with respect to the IgG for control sample which has lot more reads at the same time Knockout sample 1 has less read. So is it good to go for further analysis or i need to down sample them ?
Any suggestion is really appreciated

Usually downstream programs will take sequencing depth differences into account for you, so you only need to downsample if you get really funky results.

"if you get really funky results." i did plot using your deeptools i do see a small peak in IgG  so i feel there might be some issue.. which one would be better IgG or input im not clear .for the comparison ..any insight if IgG would be better or input ?

control is almost always better than IgG. For this reason don’t bother with IgG.

No, downsampling in general has no use as it reduces power. Peak callers such as macs2 will internally scale data to account for sequencing depth differences.
The IgG itself is only useful for peak calling as it somewhat corrects for the unspecificity of the antibody to precipitate any kind of protein rather than the protein of interest. I hear people say, and personally agree with this by experience with my own experiments, that of course IgG pulls down notably fewer DNA that the target antibody, therefore the library requires notably PCR to produce enough material for sequencing. And this then makes it questionable how representative IgG is as a control. I personally prefer chromatin input to account for the intrinsic sequence/GC bias of the underlying genomic region in the sequencing/mapping. To my knowledge there are currently few established tools to really take these controls into account for downstream analysis such as differential testing via interaction terms because the library composition is simply so different from the IPs making normalization difficult.
Your focus imho should be rather to properly normalize your samples to correct for differences is signal/noise ratio, see for example my post over at biostars on how one could do that beyond naive per-million scaling.

yes i seen that post even i wanted to use get count use deseq2 but the issue is it needs replicate and here one one of the replicate has less reads. I started with homer but after reading a bit shifted  to macs " macs2 will internally scale data to account for sequencing depth differences." then i guess no need to downsample if i understand correctly and now I will go with input..

Yes, if you use macs you do not need to do anything, the tool will correct depth internally.
One remark, the link I shared is to visualize data for e.g. profile plots, it is not the input for peak calling. For this simply use the bam files as-is.

i ran differential peak analysis using homer  getdiffretiapeak function since only one replicate from ctrl and treatment worked..now is it advisable to run with one replicate and in the  peak file argument ..so do i have merge both control and treatment peak and give or only one peak? i ran using merged both peak the result seems opposite..im not sure what is wrong..

My advise is to run a proper experiment by doing replicates. I cannot comment on Homer since I use edgeR which is more reliable as it takes replicate information into account. ChIP-seq is so noisy (by personal experience) I would never do unreplicated ChIPs.

"ChIP-seq is so noisy " im having the hard time RNA seq looks so easy now with compared to chip seq


text:
How to count the kmer occurrence in FASTA file considering overlapping and reverse complement?
I am using count_overlap() for counting the kmers from Biopython .
Does it take the reverse complement of kmer into account?
I need to count reverse complement too

Hi @Lipsa, please provide a test example and your output together with any possible error messages. This appears certainly solvable. Is this method in Biopython?

I've had a guess at what I think the question should read

Thank you for your response Michael . However , this approach is getting very very slow when FASTA is of very large size and kmers are in millions. Any help is appreciated.

There's a lot detail missing here, but at a guess could you try,
from Bio import SeqIO
from Bio.Seq import Seq
for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    reversec = record.seq.reverse_complement())
    revcount = reversec.count_overlap('AT')
    print (revcount + count_overlap('AT'))

Essentially count_overlap only counts the forward strand if this is the Biopython method. Not debugged the above.

Here's a C++ kmer counter I wrote for my own purposes. It can count canonical and reverse complement kmers from FASTA files up to length 31. Run kmer-counter --help or read the README for count options.
For longer kmers, you'd want to look at an approximation tool like Jellyfish.

I would NOT recommend using biopython for a task like this. First, it's slow, but also it's really not necessary. There are several other existing tools optimised for exactly this purpose.
My favourite is KMC. It's very fast, and by default, it reports canonical kmers (Parameter -b would turn this off). So, assuming you have KMC installed, this is how you get all the kmers from a range of occurances (1 - 10000).
# just making a directory for temprary files
mkdir tmp

# building a kmer database with k = 21,
# 16 cores and
# up to 64GB of memory.
# Counting occurances between 1 to 10000,
# `-fa` specifies that the input is fasta, you can also input other formats
# the database name will be kmer_counts
# & saving temporary data to tmp directory;
kmc -k21 -t16 -m64 -ci1 -cs10000 -fa my_seq.fasta kmer_counts tmp # run kmc

# extract kmers from the database in a text file
kmc_dump -ci1 -cx10000 kmer_counts kmer_k21.dump

I am lost here. Please help! I am using python 3. One vaex dataframe has 1million kmers and I want to check the count of these kmers in the FASTA file. How will I use KMC for this?

Yes, KMC can count kmer coverages directly from a fasta file. Just follow the code I posted there on how to use KMC to get kmers and their corresponding coverages.


text:
How do you visualize minimap2's PAF output format?
As described by the developer of minimap2:

PAF is a text format describing the approximate mapping positions between two set of sequences.

He also suggest to use either Ribbon to visualize PAF files, or to use minidot from miniasm:

minimap2 -c -x??? asm1.fa asm2.fa > out.paf
minidot out.paf > out.eps && epstopdf out.eps

What does -x??? stand for? And how do you visualize PAF files? Could you think of a way to use IGV to visualize PAF files?

As described in the minimap2 manual, -x STR is used to specify a preset option like -x map-pb for PacBio/Oxford Nanopore read to reference mapping (-Hk19).
In general, I use the -a flag to output alignments in sam format that can then be  sorted and converted to bam format, and viewed in any genome browser. The support for PAF format when it comes to visualization appears to be limited. A few options that I have found, but not tested are:

minidot
pafr
pafpy -- this one may not have the plotting capabilities though

Depending on what level of information you need, it is not difficult to just homebrew a plot in your plotting package of choice, as PAF has relevant coordinate information. Unless you are interested in some more complex features of a plotter such as IGV, I think that this is a good flexible solution.
For example, this R code makes a dot plot of the synteny between the reference and the query:
library(ggplot2)
paf = read.table("out.paf", sep="\t")
ggplot(data=paf, aes(x=V3, xend=V4, y=V8, yend=V9, color=factor(V1))) + geom_segment() + labs(x="Reference coordinate", y="Query coordinate")

The plot looks like this:

It is effectively a rather minimalist version of what pafr (as suggested by vkkodali) or minidot is doing.
(h/t Zev K.)


text:
Functional Annotation vs Functional Enrichment--which one should I use for Network Analysis?
I am currently working on a network pharmacology project where I am supposed to find out the possible diseases our drug might have an effect on. I am not a bioinformatics student, so my knowledge in the subject matter is very limited. I have all the targets and I am supposed to get the pathways associated with the targets. I used STRING for the purpose.

I got both KEGG functional annotation and enrichment for my list of targets.
I have to build a two-mode network

Question I am not very clear about

the difference between enrichment and functional annotation and
which data should be used for the next step.

If you can add more detail that might help, such as the pathway(s) your STRING analysis indicates. Your question appears vague for the decision you have and it is not clear why you don't use both for two different network analyses

The answer to your main doubt about "annotation vs. enrichment" is that you need both, but you already have one annotation.
You can only perform a functional enrichment analysis after you have a functional annotation. In your question, you state that you have KEGG annotation for your targets. That is it, knowing which target "belongs" to each KEGG module is an example of functional annotation of your targets. Several tools can annotate genes, proteins, etc. Knowing which protein is related to each drug or disease is also a functional annotation.
Now, knowing which targets are significantly important can be achieved by functional enrichment analysis.
For the network, both information can be used together. For example, you could link all elements to their functional annotation and use the available enrichment score as a quantitative parameter of importance for calculating correlations.


text:
I have used TER to break the long bonds of a chain in my PDB
I’m now not sure what I need to alter in my PDB to get it to work in leap. I know breaking the bonds turns the formerly connected residues into terminal residues.
It keeps saying that 3 of my atoms no longer have a type , how can I fix this ?

For the rest of us could you explain leap please?

Normally you have three options in a PDB:

Add a TER record
Add the  OXT atom, if your chain was cleaved or cleaved spontaneously (i.e. Asp-Pro sequence), the C-terminal would sill have a terminal oxygen. i.e. making it a terminal residue. If you have protons, you'll need to add also 3H .
Assign a different chain

If you are using forcefields as is your case, the second is the only choice. The C  atom and the N  atoms have a type that does not match because they have the wrong number of bonds —hence the type error.


text:
Moving UMI from tag into read name
I have a bam file with Unique Molecular Identifiers (UMIs) for each read present on the RX tag ( such as RX:Z:TGAGAAGGG), as expected by picard and fgbio tools.
However, several other tools (such as UMI-tools) expect the UMIs to be present at the end of the read name, such as (NB501943:73:HK333BGXC:3:12402:20313:16437_TGAGAAGGG).
Does anyone know of a tool to do this?
I have an awk solution which I'll self-answer with, but I don't know whether that is robust.

Did you check UMI-tools, this is afaik exactly what it is doing?

@ATpoint, unless I'm incorrect, UMI-tools can't do this directly. umi_tools extract moves sequences from a position in the FASTA sequence onto the read name, but in this case the reads have been added to a tag via fgbio.

My current solution is (awk is not a tool I'm very familiar with, so this may not be good awk code, there is presumably a way to get it to print the header lines unmodified.):
samtools view -H temp.sam > temp2.sam
samtools view temp.sam |
awk '{OFS = "\t"} {for(i=1;i<=NF;i++)  if ($i ~ /^RX:Z:/) {$1=$1"_"$i; gsub("RX:Z:","",$1); print}}' >> temp2.sam

This seems to work when I test it, but I don't know if I'm missing anything.


text:
Plot information stored in dataframe within a tree (ggtree)
Hello to the entire Stackoverflow community!
I'm writing to you because I'm currently building a phylogeny with ggtree and I have in parallel a table like this one:
df :
Node_number Number
50   31
1    15
13   25
35   12
40   4
5    2

and a tree such as :

and the idea is that I woule like to incorporate the df$Number information in each Node in the tree and get somethhing like I put in the picture (expected part)
tree data :
structure(list(edge = structure(c(31L, 32L, 33L, 34L, 34L, 35L,
35L, 33L, 32L, 36L, 37L, 37L, 36L, 38L, 39L, 40L, 40L, 39L, 41L,
42L, 42L, 43L, 43L, 41L, 44L, 44L, 45L, 45L, 46L, 46L, 47L, 48L,
48L, 47L, 38L, 49L, 50L, 51L, 51L, 50L, 52L, 53L, 53L, 54L, 54L,
52L, 49L, 55L, 55L, 31L, 56L, 57L, 57L, 58L, 58L, 56L, 59L, 59L,
32L, 33L, 34L, 1L, 35L, 2L, 3L, 4L, 36L, 37L, 5L, 6L, 38L, 39L,
40L, 7L, 8L, 41L, 42L, 9L, 43L, 10L, 11L, 44L, 12L, 45L, 13L,
46L, 14L, 47L, 48L, 15L, 16L, 17L, 49L, 50L, 51L, 18L, 19L, 52L,
53L, 20L, 54L, 21L, 22L, 23L, 55L, 24L, 25L, 56L, 57L, 26L, 58L,
27L, 28L, 59L, 29L, 30L), .Dim = c(58L, 2L)), tip.label = c("t5",
"t9", "t25", "t14", "t30", "t26", "t2", "t10", "t12", "t15",
"t1", "t3", "t6", "t23", "t17", "t18", "t4", "t22", "t27", "t16",
"t29", "t7", "t21", "t20", "t13", "t11", "t19", "t24", "t8",
"t28"), edge.length = c(0.912875924259424, 0.293603372760117,
0.459065726259723, 0.332394674187526, 0.65087046707049, 0.258016780717298,
0.478545248275623, 0.766310670645908, 0.0842469143681228, 0.875321330036968,
0.339072937844321, 0.839440350187942, 0.34668348915875, 0.333774930797517,
0.476351245073602, 0.892198335845023, 0.864339470630512, 0.389989543473348,
0.777320698834956, 0.960617997217923, 0.434659484773874, 0.712514678714797,
0.399994368897751, 0.325352151878178, 0.757087148027495, 0.202692255144939,
0.711121222469956, 0.121691921027377, 0.245488513959572, 0.14330437942408,
0.239629415096715, 0.0589343772735447, 0.642288258532062, 0.876269212691113,
0.778914677444845, 0.79730882588774, 0.455274453619495, 0.410084082046524,
0.810870242770761, 0.604933290276676, 0.654723928077146, 0.353197271935642,
0.270260145887733, 0.99268406117335, 0.633493264438584, 0.213208135217428,
0.129372348077595, 0.478118034312502, 0.924074469832703, 0.59876096714288,
0.976170694921166, 0.731792511884123, 0.356726912083104, 0.431473690550774,
0.148211560677737, 0.0130775754805654, 0.715566066093743, 0.103184235747904
), Nnode = 29L), class = "phylo", order = "cladewise")

does someone have an idea ?

I will have a look @chippycentra later in the week. The idea is cool. I don't really do R, but I will look into it. There are pop gen solutions to this issue but don't think that quite fits here.

You can accomplish this using the ggtree package available on Bioconductor.
First you will need to combine your tree with the data.
library(tidyverse)
library(ggtree)

ftree <- tree$edge %>%
   as_tibble() %>%
   mutate(Node_number = 1:n()) %>% # finds edge numbering
   right_join(data, by = "Node_number") %>% # find internal node associated with your node numbering
   left_join(fortify(tree, ladderize = FALSE), ., by = c("node" = "V2")) # merge tree and data

Note that this will produce a warning, but it does work.
Then plot the tree using ggtree and its geom_nodepoint layer.
ggtree(ftree, ladderize = FALSE) + # plots tree
   geom_tiplab() + # adds tip labels
   geom_nodepoint(aes(subset = !is.na(Number), size = Number, x = branch)) + # adds your data as points
   scale_size_continuous(range = c(4, 12)) # scales point sizes

You can adjust the size of the points with the scale_size_continuous function.
Data
data <- read_csv('Node_number,Number
50,31
1,15
13,25
35,12
40,4
5,2')

and tree is the same as you provided.
Ggtree also has a decent manual available at https://yulab-smu.top/treedata-book/.

thank you very much for your help :-)


text:
What datasets are available out there for prediction based on DNA sequences?
I am looking for publicly available data for a genomics deep learning project. My goal is to compare different architectures to predict biological insights from DNA sequences.
I have heard about Janggu before, but I am looking for a few datasets which are like the MNIST equivalents in genomics. Kaggle doesn't seem to have many such datasets, and I would be happy if you share your insights about potential datasets that one could use.
For example, it could be datasets that map 100-Mers to their ChIP-seq values or datasets that annotate a DNA sequence with the probability that a transcription factor would bind to it.
Choong and Lee lists a few datasets on their IEEE paper from 2017:

related:
https://www.researchgate.net/post/Are_there_any_opensource_datasets_which_contain_the_DNA_of_an_individual_and_the_diseases_asssociated_with_this_DNA

If you add more details about the type of deep learning you're evaluating we will be better equipped to direct you to resources. If you're doing something like compositional data analysis, check out this paper and the list of 45 datasets with links on the last page.

@jared_mamrot thank you. I would like to train networks that take chip-seq data and DNA sequences as inputs. The prediction layer may be any biological signal one could think of?

Genomics is a huge field; you need to be more specific. If you're looking for everything with regards to publicly-available datasets, then there's GenBank

David Reich Lab https://reich.hms.harvard.edu/datasets

Good question and the Janggu link is cool. Thanks!
Kaggle is/was heavily image recognition. There isn't an MNIST type data set in genomics, in case anyone is wondering this is a classic natural language processing (NLP) type analysis to identify handwritten single digit numbers.
There has been several FDA precision challenges, one being described here:
FDA precision challenges
This could provide something towards a 'standard'.
A key issue remains the accurate vectorisation of sequence and whether you think that is/was satisfactory in my personal opinion.
There is an old review article here

Since I have a new method to address this problem. I would like to find a dataset which is like MNIST to evaluate my idea.


text:
Why ORF13 and ORF14 of Bat Sars coronavirus Rp3 have no corrispondence in Sars2?
I'am comparing bats coronavirus with Sars2, i.e. COVID-19. I have one question...

one Bat Sars coronaviurs identify in GenBank as DQ071615.1 has the ORF13(28122..28415) and ORF14(28572..28784) proteins
Sars2 identify in GenBank as MN908947_3, has range 28284..28577
not labeled as ORF13 and the same 28734..28955 is not labeled as ORF14?

The ORF13(28122..28415) of Bat coronavirus and  PSEUDO_ORF13(28284..28577) of sars2 share 72% homology
and ORF14(28572..28784) of Bat coronavirus and  PSEUDO_ORF14(28734..28955) of sars2 share 75% homology

This seems one publication for ORF14 and ORF13 in Bat Viruses and relations in Sars2 here
The author say that ORF14 in Sars2 is all contained in N protein...
And ORF13 is know in that above link as ORF9B
and in that written is said ORF9B is implied in immunitary innate suppression...

Key to the answer is whether Rp3 is a Chinese bat virus ... I think the virus was originally part of this study https://science.sciencemag.org/content/310/5748/676 If you can tell me its country of origin .. if its Europe or Africa that would explain it.

i remember only i was surprised orf1ab was so different from sars2. And has in orf3 one 7 letter amminoacid chain that i was sarching in common with ebola (if i remember well)

The virus Rp3 stands for Rhi pea isolated in China in 2004.
The trees is here and groups with the Bulgarian sarbecoviruses and is closely related to the African sarbecoviruses. Thus although it is a Chinese sarbecovirus is not with the "Chinease originating" clade as is SARS and Covid-19.
Thats why the more mysterious ORFs ORF13 and 14 are missing, because over notable genetic distances these genes along with ORF10 are present/absent. This is a common feature of the betacoronaviruses.
In summmary this virus is notably distance from the major human pathogenic SARS/covid-19 viruses. Thus the absence of these auxilary genes isn't really surprising.
Its the very last taxa on the SARS-related viruses clade (monophyly) and is in a different group to SARS, which is the Beta-1 clade on this most excellent paper.

To answer the comment. They are outgroups, basically ancestral to COVID-19 and SARS-1 clade. So they lost these genes whilst the ingroup including COVID-109 and SARS-1 retained ORF13, 14 ... but I think lost ORF10 from memory

"Thats why the more mysterious ORFs ORF13 and 14 are missing, because over notable genetic distances these genes along with ORF10 are present/absen"... the distance of what? >70% the same

70% is a lot, post modified above.


text:
How to choose docking score cutoff?
I made a virtual screening of a large database (~5M) and want to filter ligand structures for further more accurate screening (from Schrodinger HTVS to SP resolution). My next step is to define a cutoff docking score, but it seems arbitrary to choose one. Are there any established methods for filtering the ligands out?

It is utterly arbitrary and depends on your resources.
We can all agree your cutoff would be a lot lower during the ordering step than at followup in silico analyses, but even then if you were ordering the top N at \$100-250 each, should you spend \$1,000 or \$10,000?
In the in silico step you have real time and CPU-time to factor in making this call even blurrier.
But there are some things to put those scores into context:

The units are kcal/mol, i.e. these are predicted ∆∆G of binding. Here is an explanation of what does a kcal/mol mean for docking — 1 kcal/mol is a hydrogen bond or the average collision energy of water at 37°C.
If you have a crystal structure with a known ligand, re-docking and cross-docking that would give you a benchmark of what you should expect ——keeping my mind that natural ligand are horrid and don't obey Lipinski's rule. The natural ligand will wander and to assess how good the docking went RMSD is used, but for this I mean looking at its GlideScore as if it were any other compound.
Rescore with a different tool —Gold etc.

Ligand efficiency
An important thing to remember is that ∆∆G is size dependent —a huge peptidomimetic will score more than a smaller molecule that hits the hotspots perfectly. Therefore, ligand efficiency (LE), ∆∆G over heavy atoms (reviewed in detail here) is a good metric to prevent huge molecules dominating spuriously. If your library has a variety of different sizes, it is best to sort by LE than ∆∆G alone. In fact, if you get a few small hits you can merge/link them or expand on them, whereas if you get nothing because you gambled on larger hits you have to start from square one. I would argue that this is true for repurposing screens and the low success of repurposing VLS despite the enthusiasm is due to the neglect of this.

Great answer. Just to add context to the kcal/mol point; these will only be close to accurate if the parameters and assumptions perfectly mimic the chemistry, and even then the biochemistry can complicate things further as the linked article points out.

@James Your link did not stick. But thanks —you reminded me that I did not mention ligand efficiency!!

Sorry Matteo, I was meaning to refer to the link in your question explaining kcal/mol in docking.


text:
Convert sequence in MS doc to fasta or genbank
How can I convert a sequence provided in a Microsoft doc file into a fasta or genbank format?

The CrazyDoc Python package can convert Genbank/Fasta/Snapgene/MSDoc sequences into Biopython records, which can be saved as Genbank or Fasta.
A web interface / demo is provided at EGF CUBA (Collection
of Useful Biological Apps): Convert Sequence Files.
Disclaimer: I'm the current maintainer of CrazyDoc.

A fasta file is just in plain text format. Pandoc can convert .docx format to plain text format:
pandoc -i input.docx -t plain >output.fasta


text:
How to fix NameError?
I am working on a project and am having issues with the following code that I have written in nano:
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
       print (">"+list.group(1))

This is the error I receive when I execute my command on command-line:
File "mpo.py", line 7, in <module>
    gene_id = myfile.id
NameError: name 'myfile' is not defined

I have a fasta file with the format
>H149xcV_Fge342_r3_h2_d1 len=210 path=[0:0-206]
ACTATACATGAGGAGAACATAGAACAAAAATGGGACCATAGATATATAACAATAGAAGATATAGAGAACACAATAGACAACTTATTAGGAAAGAGGTGTGTCGTCATGGAGCTGATGTTCGAGGATACTTTGCATGGTCATTCTTGGATAATTTTGAGTGGGCTATGGGATACACCAAGAGGTTTGGCATTGTTTATGTTGATTATAAGAACGGGC

>H149xcV_ytR1oP_r3_h2_d1 len=306 path=[0:0-207]
ATTAGAGTCTGAGAGAGTCTTGATTTGTCGTCGTCGAGAAATATAGGAGATCTGATTAGAGGAGAGAGCGGCCTAGGCGATGCGCGATATAGCGCTATATAGGCCTAGAGGAGAGTCTCTCTCTTTTAGAAGAGATAATATATATATATATATGGCTCTCCGGCGGGGCCGCGCGAGAGCTCGATCGATCGATATTAGCTGTACGATGCTAGCTAGCTTATATTCGATCGATTATAGCTTAGATCTCTCTCTAAAGGTCGATATCGCTTATGCGCGCGTATATCG

I wish to reformat my file so that it only provides me with the unique gene id's for genes with a length greater than 250 bp.
My desired output is:
>H149xcV_Fge342_r3_h2
>H149xcV_ytR1oP_r3_h2
>H149xcV_DPN78333_r3_h2
>H149xcV_AgV472_r3_h2
>H149xcV_DNP733_r3_h2

note that after your edit, your code now does not match your error. You should update the error to the new error, otherwise we can't help you very well.

Hi @MaximilianPress I have updated the error!

I updated the post to attempt to trap the bugs.

You might consider changing the code to make the whitespace a little more consistent (as it was, it failed when i tried to run it).
Also, why are you using regular expressions? You have already extracted the gene_id using Biopython. The following works for me:
(updated to handle the underscore delimiting)
from Bio import SeqIO
import sys
import re

fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile.seq) > 250:  # note .seq here!
        gene_id = myfile.id
        # unclear why using regex
        #match_list = re.match('H149xcV\_\w+\_\w+\_\w+',gene_id)
        #print (">"+match_list.group(0))
        print("_".join(gene_id.split("_")[0:-1]))

Hi @MaximilianPress , how would I change the regex so that it only provides me with >H149xcV_Fge342_r3_h2 and not>H149xcV_Fge342_r3_h2_d1. I would not like the very end part that consists of d1. So basically, how can I format it so it it only provides this >H149xcV_Fge342_r3_h2. If I am using regex would I not need to use gene.id?

@andnowmywatchbegins wants to remove the stuff after the space in the fasta header. In phylogeny this is either removed or spaces substituted because the space in phylip format denotes the start of the sequence

@Michael I've tried reformatting the entire regular expression and I am still getting the full >H149xcV_Fge342_r3_h2_d1, is there a way to only receive >H149xcV_Fge342_r3_h2 as in output? Thanks!

Can you put this as a separate question please? .. Yes you will but thats a different part of regex

@Michael stuff after the whitespace is already stripped out by biopython when you use the .id. this is one of the more annoying behaviors of biopython IMO, but in this instance it is useful. I've confirmed that my example code does this. It admittedly does not handle the underscore delimitation, which I missed. But that can be done with split much more simply.

@andnowmywatchbegins why not just split on _ and only keep as many _-delimited fields as you want? Regex will work for this in theory but it's annoying to debug in cases like this.

@MaximilianPress I've put a response here .. https://bioinformatics.stackexchange.com/questions/14913/how-to-obtain-desired-output

@MaximilianPress "if len(myfile.seq) > 250:  # note .seq here!" -- .seq is actually not needed there

@Chris_Rands ah, good to know... another confusing feature of biopython!

You wrote,
mylist = re.match(H149xcV\_\w+\_\w+\_\w+)

try...
myregex = re.compile(r'H149xcV\_\w+\_\w+\_\w+')
mylist = myregex.match(myfile.id)

or
mylist = re.match(r'H149xcV\_\w+\_\w+\_\w+',myfile.id)

Watch your indentation it should be 4 spaces. Note you are losing the last '_\w+', but don't think it matters.
I think it is better to use,
mylist, _, _ = re.split(r' ', myfile.id)

The problem with your code is list will just return true ... and that is not what you want. If you are using split carefully watch the number of ' ' in your code you can also use r'\s'
You can also set maxsplit as a parameter, which might be useful.

myfile not defined.
Please try,
fasta_file = (sys.argv[0])

Generally your argv command isn't cool and not just [0] versus [1]. I forget the argv module everyone uses .. there's two of them, they are a bit fiddle but are very cool of options. You way in is okay, okayish, sort of okay ... Anyway you appear to be requesting the second input into the commandline, wherea you are probably placing the fasta file in the first position.
Also
print (fasta_file)

Check something happens that you picked up the file.
The other way in is,
mypath = '/User/username/location/fasta.fa'
for myfile in SeqIO.parse(mypath, "fasta"):

where mypath is a OS X style, you're probably using Linux so its... /home/username/pathtofile/fastafile.fa (something like that) or
for myfile in SeqIO.parse('/User/username/location/fasta.fa', "fasta"):

I forgot .... don't ever use list as a variable, by mylist is fine. Again I've not debugged my code.

Hi @Michael I tried doing that, but received a NameError indicating that the name 'myfile' is not defined. I appreciate your help!

@andnowmywatchbegins can you update your answer with this information?

list is used in Python for creating lists. I would recommend naming this variable something else.

list .... Yes!! Use mylist as a variable. instead


text:
How to obtain desired output?
I am working on a project using the following command within nano:
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match(r"H149xcV_\w+_\w+_\w+", gene_id)
       print (">"+mylist.group(0))   # edited from list to mylist

and its providing with the following outout:
>H149xcV_Fge342_r3_h2_d1
>H149xcV_bTr423_r3_h2_d1
>H149xcV_kN893_r3_h2_d1
>H149xcV_DNp021_r3_h2_d1
>H149xcV_JEP3324_r3_h2_d1
>H149xcV_SRt424234_r3_h2_d1

How can I change my command so that it provides me with the desired format and only the UNIQUE gene id's:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

its easy replace match with sub but please stop using list as a variable ... mylist is fine.
This might work
   mylist = re.sub(r'H149xcV_\w+_\w+_\w+', gene_id)

otherwise
myregex = re.compile('_\w+\s+.*')
fastaid = myregex.sub('', myfile)

OR from @MaximilianPress
myregex2 = re.compile('_\w+\n') # or myregex2 = re.compile('_\w+$')
fastaid2 = myregex.sub('\n', myfile) # or fasaid2 = myregex.sub('', myfile)

The above will work ... As with all my code I never proof it ..

Hi @andnowmywatchbegins thats easy just use dictionary to remove the duplicates ... its been answered several times on the site. Please remember accept/upvote stuff.

from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_\w+_\w+_\w+", gene_id)
            print (">"+mylist.group(0))

let me know if it works I am learning too


text:
How to remove duplicates from a fasta file using python
I am using the following command:
from Bio import SeqIO
import sys
import re

fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
if len(myfile.seq) > 250:
    gene_id = myfile.id
    mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
    print (">"+mylist.group(1))

And this is providing me with duplicates of the same gene:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2
>H149xcV_Fge342_r3_h2
>H149xcV_Fge342_r3_h2

How can I reformat my command so that I only receive unique gene id's:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

Why do you need to use python for this?

I am using nano within commandline to write my python command. @RamRS

I'm asking why, not how. There are command line tools that de-dup fasta, so why reinvent the wheel?

@RamRS The only command-line tools that I am aware of are sort|uniq but I am not sure as to how I would use it? I must apologize, I am kind of new to bioinformatics, so I am still getting the hang of it. I know I can output my python into another file, but is there a way to obtain the unique id's afterwards? Thanks!

You can use  seqkit rmdup instead reinvent the wheel. https://bioinf.shenwei.me/seqkit/usage/#rmdup

@zorbax I think that would be a good answer...

I wouldn't use Python for this, myself. Instead you can use the FastaToTbl and TblToFasta scripts I have posted before, and pipe to standard *nix utilities:
FastaToTbl file.fa | sort -u | TblToFasta > file.uniq.fa

Or, if you don't have GNU sort:
FastaToTbl file.fa | sort | uniq | TblToFasta > file.uniq.fa

Alternatively, you can also do something like this:
awk -v RS='>' -v ORS='>' '++a[$0]==1' file.fa | sed 's/>$//'

For example, given this input file:
$ cat file.fa
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA

You can do:
$ awk -v RS='>' -v ORS='>' '++a[$0]==1' file.fa | sed 's/>$//'
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT

Both of the above approaches have the advantage of not needing to store the entire file in memory.

There is also an application at UCSC utilities that converts Fasta to Tab separated file. It's called faToTab and can be downloaded from http://hgdownload.soe.ucsc.edu/admin/exe/

from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
            print (">"+mylist.group(1))

Yeah, I would do something along these lines too. Though, I would use a set or dict insead of list if the fasta file is big...

Perhaps itertools might be an option? From the itertools page / moreitertools project
https://docs.python.org/3/library/itertools.html#itertools-recipes
def unique_everseen(iterable, key=None):
    "List unique elements, preserving order. Remember all elements ever seen."
    # unique_everseen('AAAABBBCCDAABBB') --> A B C D
    # unique_everseen('ABBCcAD', str.lower) --> A B C D
    seen = set()
    seen_add = seen.add
    if key is None:
        for element in filterfalse(seen.__contains__, iterable):
            seen_add(element)
            yield element
    else:
        for element in iterable:
            k = key(element)
            if k not in seen:
                seen_add(k)
                yield element

I would use a dict so that the lookup time is constant. If you are creating an array with all the unique gene IDs, then checking if a gene is in that array, your runtime is O(N2). When you create an entry in the dict, the time to access is constant, and if the entry is already present, it will just overwrite it. Something like this should do it for you:
from Bio import SeqIO

genes = {}
fasta_file = sys.argv[1]
for entry in SeqIO.parse( fasta_file, "fasta" ):
    genes[ entry.id ] = entry.seq

Then you can write a new, unique fasta file:
with open('unique_genes.fasta', 'w') as out:
    for key,value in genes.items():
        out.write( '>%s\n%s\n'%(key,value) )

Your question isn't well formatted, for example list is used as a variable name where it should be reserved as a keyword...
print (">"+list.group(1))

However, the answer you're after is going to be something like:
print (">"+set(name_of_list.group(1)))

Anytime you want to get the unique elements of a list, just wrap it in set()

I use the following code. It doesn't only remove repeated sequences but similar sequences (with a quite simplistic approach).
I use it just sporadically and I guess it can be largely depurated. Remove the debug prints if these annoy you.
import sys
import requests
import shlex, subprocess
import random
from Bio import SeqIO
import editdistance

if len(sys.argv)<3:
    print("Uso: "+sys.argv[0]+" input_fasta output_fasta")
    sys.exit(0)

savefile = sys.argv[2]
content_file = sys.argv[1]

class Sequence:
    id=""
    sequence=""

def are_similar(seq1, seq2):
    l1 = len(seq1.sequence)
    l2 = len(seq2.sequence)
    if abs(l1-l2) < l1/40:
        d = editdistance.eval(seq1.sequence[0:250], seq2.sequence[0:250])
        #d = editdistance.eval(seq1.sequence, seq2.sequence)
        print(seq1.id+" "+seq2.id+" -> "+str(d))
        #if d < l1/20:
        if d < 50:
            return True
    return False

def find_repeats(sequences, seq, repeated):
    b=False
    for seq2 in sequences:
        if(b and not seq2.id in repeated):
            if are_similar(seq, seq2):
                repeated.add(seq2.id)
        if seq2.id==seq.id:
            b=True

sequences=[]

res=""

for record in SeqIO.parse(content_file, "fasta"):
    seq=Sequence()
    seq.id=record.id
    seq.sequence=record.seq
    sequences.append(seq)

repeated_sequences=set()
for seq in sequences:
    find_repeats(sequences, seq, repeated_sequences)

print(repeated_sequences)

def mySortFunc(e):
    return len(e.sequence)

sequences.sort(key=mySortFunc)

for seq in sequences:
    b= seq.id in repeated_sequences
    if not b:
        res+=">"+seq.id+"\n"
        res+=str(seq.sequence)+"\n"

output_file = open(savefile, "w")

output_file.write(res)

output_file.close()

EDIT: I see that's not exactly what OP is asking ...


text:
Any suggestions for cultivar identification using SSR (simple sequence repeat) markers
This is the data I have now: 30 simple sequence repeat (SSR) markers for 80 cultivars of cucumber. 10 of the 80 cultivars belong to one cultivar (let's say A).
My goal is to classify an unknown cultivar into A or "not A" using tge 30 SSR markers.
I'm think of considering it as a classification problem (A vs non-A) and use machine learning method to build a model using the SSR markers as features. But the problem is that A cultivars don't have enough number of samples.
Do you have any suggestions which statistical method(s) I can try to solve this problem? Thanks in advance.

Use Fst (migration statistic) such as the Fstat program to generate a distance matrix and solve it via clustering, such as UPGMA or neighbor-joining. This analysis is particularly useful for diploids. I presume you have heterozygotes in your analysis. You can switch Fis (inbreeding) if you have excessive homozygosity.
The other way in is a haplotype network. When I get a minute I'll post the link to the programs. Machine learning isn't needed in this instance, although you could do a PCA style analysis, i.e. unsupervised learning. Unsupervised learning doesn't need a large sample size.

The other solution is a stepwise mutation phylogeny, but this is for microsatellites (one to four repeats). This is called Delta Mu squared and describes the method of repeat unit mutation. It is a distance method.
Again when I have a moment I will supply weblinks.

Not to discount @Michael's answer (unsupervised learning is very handy for descriptive analysis of this kind of problem), but classification based on microsatellite data should not be too difficult. For papers with examples, you can see e.g. this and this. They seem to apply various algorithms, you might be able to learn something from those papers.
More generally, I bet that you could get pretty far with just (multinomial) logistic regression for your outcome variables.
Of course, this depends to a certain extent on how you code the data and what the input data look like. Microsatellite allele coding and evolutionary/error models are a field in themselves.


text:
Reading nested Map data structures in WDL
I'm having trouble getting WDL's read_json to work as expected when trying to create a Map[String, Map[String, String]] data structure. If I explicitly declare the data structure it works, as shown below (the gm variable):
workflow TestWorkflow {

    String genome_choice
    String output_name = "xyz.txt"
    Map[String, Map[String,String]] gm = {
        "A": {"gtf":"A.gtf", "idx":"A.idx"},
        "B": {"gtf":"B.gtf", "idx":"B.idx"}
    }

    call something {
        input:
            genome_file_map = gm[genome_choice],
            outname = output_name
    }

    output {
        File fout = something.fout
    }
}

task something {

    String outname
    Map[String, String] genome_file_map

    String idx = genome_file_map["idx"]
    String gtf = genome_file_map["gtf"]

    command {
        echo ${idx} >> ${outname}
        echo ${gtf} >> ${outname}
    }
    output {
        File fout = "${outname}"
    }
}

However, instead of explicitly writing gm, I would like to read that from a file. If I then create a JSON-format file:
{
  "A": {"gtf":"A.gtf", "idx":"A.idx"},
  "B": {"gtf":"B.gtf", "idx":"B.idx"}
}

and submit that as an input file
workflow TestWorkflow {
    ...
    File resource_map_jsonfile
    ...
    Map[String, Map[String,String]] gm = read_json(resource_map_jsonfile)

it fails to parse, giving me:
ERROR: gm is declared as a Map[String, Map[String, String]] but the expression evaluates to a Object:

    Map[String, Map[String,String]] gm = read_json(resource_map_jsonfile)

If I change the type declaration to Object, then all the bracket-style indexing fails. So that's out.
Curiously, if I instead declare gm as a Map[String,String] it then fails giving:
java.lang.IllegalArgumentException: No coercion defined from wom value(s) '{"gtf": "B.gtf", "idx": "B.idx"}' of type 'Map[String, String]' to 'String'.

so clearly it is parsing the JSON file correctly, as it identifies that nested Map[String, String]. I can probably find ways around this (using other read_* functions), but that's less than ideal.


text:
How to sort multiple FASTA files based on their content?
I have around 10,000 FASTA files of Influenza A virus.
These files contains sequences of each of the 8 segments of the viral genome and I want to separate each genome segment into a different folder and subdivide each directory into different files based on their geographic region of isolation.
In each FASTA file for each segment,the first line has the segment number, for example.

KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

To be clear I want

file1.txt has segment 1
file2.txt has segment 2
etc ...

I want to ensure all segment 1 sequences are grouped into one folder, and each file is grouped according to its geographic origin. The geographic groupings are mirrored for all 8 segements and each placed into their own directory.

Cross-posted at https://stackoverflow.com/questions/64586925/how-to-sort-multiple-fasta-files-based-on-their-content?noredirect=1 (closed now).

The general idea is:

iterate over all files

read the first line of each file, extract the segment bit
move the file based on that

In Python, this can be done in many ways. Here’s one way:
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “segment X” in their FASTA ID, will trigger an error, as will the usual IO failures (missing permissions, etc.).

But most people would probably use a shell scripting language instead of Python to solve this. For example, here’s something equivalent in Bash (untested):
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ -n $segment ]]; then
    mv "$file" "segment-$segment"
    done
done

Very cool code.

One "out there" approach you could take is to create a fasta file containing your regions of interest (as separate sequences), and then map your fasta files to that reference fasta file, generating BAM/SAM output. After that, position-sort the BAM file, and you've now got a content-sorted list of sequences.
The process to do this would be something like this:
minimap2 -a regions.fasta input1.fasta input2.fasta ... | samtools sort > mapped.bam
samtools fastq mapped.bam > sorted.fasta


text:
Reading fasta sequence regions with faidx
Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.
To do it in c++ i use the library htslib and the command faidx_fetch_seq() which is similar to read specific regions using samtools faidx.
The problem is when im reading in the reference genomes, there are of course multiple regions consisting only of NNNNN (and so on). Since i only need the actual DNA sequence i wish to read in the entire file but removing the regions with NNNN.
I have tried reading in the regions in c++ and then just creating an if-statement to check the actual bases consist of N. This works but is incredibly slow.
So my question is: do anyone know how to only read in the actual DNA sequences, removing the NNN regions? either with samtools or with the c++ funciton faidx_fetch_seq().

Can you put "incredibly slow" in context? Slow is different for everyone, are we talking seconds, minutes, hours? Can you post some of your code that you wrote for this?

I just tested with a tiny program on a 400Mbp genome, reading in all the sequence while at the same time removing Ns. Took ~10.8 seconds, is that too slow? There are certainly ways to optimize this but you will be mostly bound by disk speed at this point, checking for Ns in a string of characters is not your bottleneck

Perhaps they are working with plant genomes, some of which will be unlikely to fit into the system memory of all but the most decked-out systems.

If you can't read in the entire sequence for a given reference genome:
(1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find Ns in your FASTA file). (2) Use bedops --complement with this mask file and the output of UCSC fetchChromSizes to generate a BED file containing regions that are unmasked. (3) Repeat for other reference genomes, as needed.
You can now loop over the unmasked regions in this file, making calls to samtools faidx or faidx_fetch_seq(). These calls will return unmasked sequence.
For performance or convenience, you could read chunks of sequence into a std::deque or std::vector buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.


text:
Adding multiple dataframes to compare against genotypes - DATE/LOCATION - geotemporal incidences of SARS2 genotypes
I'm trying to reuse some code I have previously used to add extra parameters.
In the initial code, it was used to give a 1/0 counter if a specific mutation was present (first code below)
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')

for lineage in lineages:

    df = pd.read_csv('Lineages/' + lineage).dropna()
    df.columns = ["lineage"]
    df['lineage'] = [x.split('|')[0].split('_')[-1]
                     for x in df['lineage']]

    for mutation in mutations:

        mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
        mutation_header = mutation.replace(".txt", "")
        mutation_df.columns = [mutation_header]
        mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                          for x in mutation_df[mutation_header]]

        mutation_set = set(mutation_df[mutation_header])

        df[mutation_header] = [1 if x in mutation_set else 0
                               for x in df['lineage']]

    df.to_csv('Results9/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

However, I now wish to add extra parameters such as date of sample and location of sample.
The first code produces this type of counter system:
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

This code is fine and working, but if I try and add the second parameter (date/location) with the following code:
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:/Users/jones/Downloads/SARSCOV2/"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')
dates = os.listdir('Dates')
locations = os.listdir('Locations')

for date in dates:

    date_df = pd.read_csv('Dates/' + date).dropna()
    date_df.columns = ["Date"]
    date_df['Date'] = [x.split('|')[-1].split('/')[-1]
                       for x in date_df['Date']]

    for location in locations:

        location_df = pd.read_csv('Locations/' + location).dropna()
        location_df.columns = ["Location"]
        location_df['Location'] = [x.split('/')[1].split('_')[0]
                                   for x in location_df['Location']]

        for lineage in lineages:

            df = pd.read_csv('Lineages/' + lineage).dropna()
            df.columns = ["lineage"]
            df['lineage'] = [x.split('_')[-1].split('|')[0]
                             for x in df['lineage']]

            for mutation in mutations:

                mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
                mutation_header = mutation.replace(".txt", "")
                mutation_df.columns = [mutation_header]
                mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                                                for x in mutation_df[mutation_header]]

                mutation_set = set(mutation_df[mutation_header])

                df[mutation_header] = [1 if x in mutation_set else 0
                                       for x in df['lineage']]

            df.to_csv('Results_3/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

I retrieve this data, where the counting is maintaibned, but the new parameters (DATE/LOCATION) are not added
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

In short, my end goal is to:

Add Lineage ID|Location of Sample|DD/MM/YYY|Mutations| -- as the headers
still maintain the mutation counter
Use the mutation counter as a guide to determine where and when the genotypes (mutation combinations) arose.

The files where I am retrieving the parameters looks like this also:

hCoV-19/Singapore/4/2020|EPI_ISL_410535|2020-02-03

hCoV-19/USA/WA13-UW9/2020|EPI_ISL_413601|2020-03-02

hCoV-19/USA/WA-UW142/2020|EPI_ISL_416680|2020-03-11

hCoV-19/USA/WA-UW143/2020|EPI_ISL_416681|2020-03-11

hCoV-19/USA/WA-UW144/2020|EPI_ISL_416682|2020-03-09

hCoV-19/USA/WA-UW145/2020|EPI_ISL_416683|2020-03-15

hCoV-19/USA/WA-UW146/2020|EPI_ISL_416684|2020-03-14

hCoV-19/USA/WA-UW147/2020|EPI_ISL_416685|2020-03-15

hCoV-19/USA/WA-UW148/2020|EPI_ISL_416686|2020-03-14

A very laborious, and poorly explained question; I apologise in advance!

I have realised an error in my coding - however, this only allows the initiation of the mutation counter. I am still unable to add the extra dataframes. Please see the editted code.

Would you be able to share two or three of your input files?

https://docs.google.com/document/d/1fV6wiYKkTngCFOeOwcfl1JGbhpY_VpWrSpGsdTfy-bw/edit?usp=sharing
Lineage A
https://docs.google.com/document/d/1qCNtzRI6EV2GNn3BPuP6EtuB5VWkFqgmGuLT4z_XqoY/edit?usp=sharing
A222V Mutation

Sorry @BastianSchiffthaler I realise those two files may seem lacking. I was in a rush for work while sending these. So the ‘Dates’ and ‘Location’ input files are identical to the ‘Lineage’ files. I merely copied these, as the code didn’t like trying to make multiple separations from the same folder for some reason

Just to clarify, when I asked you to share some of your files, I meant to share them in a way that I can reproduce your issue from start to finish. You read from at least four different files in a quadruple for-loop and I'm pretty sure your goal could be achieved much more straightforwardly than your current implementation.
Without a reproducible example, I can offer a best-guess on what's happening: you create date_df and location_df, but after you instantiate them, you never do anything with the objects. Your final write happens from df, so you'll need to add those columns, e.g.:
df["Date"] = date_df["Date"][0]
df["Location"] = location_df["Location"][0]

The above is just an example, you'll need to make sure the correct value gets assigned to your df object instead of simply the first one.

Again, apologies I was late running for work and rather skimmed the post. So would this be implemented before the mutation_set loop or between it? Or before the df[mutation_header]? I assume the latter

I assumed the latter and it’s worked perfectly? Thanks again!


text:
How to sort multiple FASTA files based on their content?
So I have around 10,000 FASTA files of Influenza A virus. These files contains sequences of each of the 8 segments of the viral genome and I want to separate these files into different locations based on the content of these FASTA files.
In each FASTA file, first line has the segment number.
Eg.
>KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

Okay imagine you have 10,000 text files in one folder and there are 8 types of these text files based on the content of them. Eg. file1.txt has segment 1 written in it, file2.txt has segment 2 written in it upto segment 8. So I want to group the files containing segment 1 together in one folder, segment 2 in other folder together and so on.

I edited your question to make it more readable. Check if I interpreted it correctly and change if not.

Thank you, It's my first time using stackoverflow, so sorry about that.

My genetics knowledge is quite basic. Which is the part you want to sort by? (I assume you have several folders with that name where the corresponding files should go?)

We can treat these fasta files as basic text files containing sequences. So i have 10,000 text files and i want to sort them into different folders based on the segment sequence these files have. Eg. One text file will have segment 1 sequence only, and so on.  In total there are 8 segments, so 8 folders.
I hope i am clear now.

What's the "segment sequence" part? KM368312.1?

This is accession number, but each text file will have a unique accession number.
But i want to group together all the segment 1 files together and so on.

I won't guess again. What exactly is the wanted part?

Okay imagine you have 10,000 text files in one folder and there are 8 types of these text files based on the content of them. Eg. file1.txt has segment 1 written in it, file2.txt has segment 2 written in it upto segment 8. So I want to group the files containing segment 1 together in one folder, segment 2 in other folder together and so on.

@KonradRudolph: why? It's about textfiles and sorting them. That the content of the files happens to be about genetic, doesn't change that (and is completely irrelevant to the question).

@Stephan There’s naturally overlap between the two sites but there are dedicated tool and libraries for working with FASTA files. But, furthermore, the current question is completely off topic on Stack Overflow because it is a “please write code for me” question. On bioinformatics.SE we’re generally more tolerant of such questions for specific bioinformatics-related tasks.

The close notice is confusing because each close voter found a different reason to close this - you are vaguely violating several of the guidelines for What topics can I ask about here. Probably review the [help] and in particular also How to ask. If you can [edit] this into a well-scoped question about a Python program you have written, with details about what you have already tried and how it failed, it can certainly be salvaged and reopened.

That's quite easy with cmd:
@echo off
for /l %%i in (1,1,8) do (
  for /f "delims=" %%f in ('findstr  /mc:" segment %%i " *.txt') do (
    ECHO move "%%f" "segment %%i\"
  )
)

findstr /mc:"string" *.txt searches each matching file for the string and outputs the filename, if the string was found. The for /f loop around catches those filenames, so you are able to move them. The for /l loop just iterates from 1 to 8.
See for /? and findstr /? for more information about those commands.
You probably want to adapt the filemask *.txt to your needs.
Note: I "disarmed" the move command for troubleshooting. Remove the ECHO when you are sure it does what you want.

I know this is tagged cmd but it’s very unlikely that OP is actually using Windows when working with FASTA files.

The general idea is:

iterate over all files

read the first line of each file, extract the segment bit
move the file based on that

In Python, this can be done in many ways. Here’s one way:
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “segment X” in their FASTA ID, will trigger an error, as will usual IO failures (missing permissions, etc.).

But many people wouldn’t use Python to solve this but a shell scripting language. For example, here’s something equivalent in Bash (untested):
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ "$segment" != '' ]]; then
        mv "$file" "segment-$segment"
    done
done


text:
How to download sequence alignments against UniProtKB sequence database for a particular Pfam Id using R?
I am trying to download sequence alignments for families against the UniProtKB sequence database as shown in the example below using R code. Can someone kindly suggest me how to?
https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1

The following code served my purpose.
library(readr)
pfam_url <-
  'https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1'
seqs <-
  read_table(
    pfam_url,
    col_names = FALSE,
    col_types = cols(X1 = col_character(), X2 = col_character())
  )