text:
Too many requests error when querying NCBI
I keep getting a "too many requests error" when querying the NCBI SRA database, even though I'm running less than 10 requests per second, and I have an API key, which supposedly should allow me to run 10 per second.
Here is my (Python) code:
import subprocess
import concurrent.futures
import time

bioprojects = [
    "PRJNA644722",
    "PRJNA644892"
]

all_metadata_generators = []

def fetch_metadata(bioproject):
    return subprocess.run(f"esearch -db sra -query '{bioproject}[bioproject]' | efetch -format runinfo", shell=True, capture_output=True, text=True)

with concurrent.futures.ThreadPoolExecutor() as executor:
    for i in range(0, len(bioprojects), 10):
        all_metadata_generators.append(executor.map(fetch_metadata, bioprojects[i:i+10]))
        time.sleep(1)

for metadata_generator in all_metadata_generators:
    for metadata in metadata_generator:
        print(metadata)

Ahh... this was often rumoured and there whispers on the wind of departmental bans - just rumours. You bang on with the sleep, I would just increase it until the signal detection police are happy. If you use a try: exception loop you could even try and attempt to attenuate the sleep time to just about the detection limit. Goodluck

@Michael thing is here sleep is technically not even doing anything since I only have 2 bioproject ids. So I can't even do 2 requests in more than one second without problems

where do you have the API key? In your .bashrc file? I would run the same command from the command line to check if you still get the 429 error. If you don't then you may have to specify the API key in the python command again.

@vkkodali how do you specify it on the command line? I couldn't find it in the e-utilities documentation

@vkkodali ah I see now that actually the API key was not set. I set it from Biopython (and I'm not using Biopython in my code) so obviously it wasn't going to be used above. Silly me!

@vkkodali but how do you pass it to the esearch / efetch commands?

As @vkkodali pointed out, the problem turned out to be that you shouldn't use multithreading.

The execution of your EDirect commands may be limited unless you are using an API Key, which you can create from your MyNCBI account. Your API Key can be used in one of several different ways.
In an argument to an EDirect command:
esearch -db nuccore -query 'some query' -api_key 12345

In a UNIX environment:
You can either run this just before you are running some edirect commands or add the following line to your .bashrc file.
export NCBI_API_KEY=12345

In the environment, from a Perl script:
$ENV{NCBI_API_KEY} = "12345"

did as you said, yet I still get "too many requests" for quite a few

Even went as far as setting the sleep time to 10, and the step size for the for loop to 1, and still get too many requests

I typically don't use Python for this sort of thing. I just run the command in a bash shell. I suspect the reason for the too many requests error is because you are using multithreading. Does your code simultaneously run 10 instances of the fetch_metadata function? May be try without the multithreading to see if it works?

yes that worked. Interesting, thanks!


text:
Extract only the amino acid residues on the surface of protein structure data
I'd like to use protein structure data from PDBs visualized by Chimera to do two things.
I'd be happy to learn the appropriate tools and methods.
・Mapping information about amino acids present on protein surfaces to sequences
・Find the charges of exposed amino acids on the surface and apply color gradations to the amino acid molecules according to their values
The second one in particular looks like it could be done with Chimera's "Coulomb surface coloring". But I'm not sure what kind of value to enter underneath the color and how to find it, which is needed to add a gradient.
Sorry for my bad English ...

Hi @farinelli what is sospecial about Chimera? For example, PyMol can be used following a surface charge calcuation

Finding solvent accessibility of residues is straightforward, however, the charge of a residue is trickier. An atom has a formal charge and a partial charge (Gasteiger-Marsili charge). Protons are not seem in protein crystals, so which protonation state a surface residue has needs to be calculated. The Advanced Poisson-Boltzman Solver for example is implemented in several tools, such as PyMOL or pdb2pqr.

source: Stack Exchange

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

sort: votes


text:
DiffBind, diffferentially binding site
I have data for 3 histone marks (2 for silencing and 1 for activation) each mark has three replicates.
when I run the diffBind package I have three contrast:
Silence1-Activation
Silence1-Silence2
Silence2-Activation

When I run dba.plotBox for any of contrasts I will have the same distribution of reads.like below:

I even tried to start with 1000 peaks of each one that showed higher fold change but again got the same pattern for plotbox. How can I see which peaks belong to which histone marks?

most likely the data frame is wrangled wrongly. Even the outliers are the same. This is almost impossible

I got the answer from Rory Stark. the samples should have a unique ID. (https://support.bioconductor.org/p/132502/#132514)

Thank you for your feedback @Mariam

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


text:
Display a score or graph to visualize the degree of conservation of each residue from alignment data
I would like to display a score or graph to visualize the degree of conservation of each residue from the amino acid alignment data.
If possible, I'd like to extract the parts of the game that scored above a certain level.
Is there a tool or source code (python) to do that?
Sorry for my bad English ...

Why not use one of the many pair-wise alignment web tools? They tend to have nice graphical outputs like that.

Thank you, which one is the best?

Just choose one, there is no best.

OK, I'll use something.

I would suggest a logo plot, e.g.:

I would take @Devon Ryan's suggestion to create a multiple sequence alignment (MSA), and plug that into that logo plot tool (or others, there are many like it).

Thank you, Max!! How do you think jaview[https://www.jalview.org/About]?

I have never used JalView but the tutorial says that it has this feature (http://www.jalview.org/tutorial/TheJalviewTutorial.pdf, p.44). I like to use smaller, simpler tools but if you like JalView it should work for making a logo plot.


text:
What is the pvalue from anova(glm1, glm2, test = “Chisq”) telling me in R when the glms have a negative binomial distribution?
Having constructed two nested negbin glms in R, one has a dropped factor, We can do an anova with test = "Chisq" to see if there is a significant difference between the models. But what does the Pvalue as a result of this actually tell me?
anova(mnegbin1, mnegbin2, test = "Chisq")
Likelihood ratio tests of Negative Binomial Models

Response: count
                        Model    theta Resid. df    2 x log-lik.   Test    df      LR stat. Pr(Chi)
1 origin + substrate + sample 18.46502      1232       -2459.886
2          substrate + sample 18.46502      1232       -2459.886 1 vs 2     0 -1.500666e-11       1

Hi @Lamma, welcome to bioinformatics. I commented on a similar question on cross-validated. https://stats.stackexchange.com/questions/476492/is-is-better-to-compare-the-log-likelihoods-between-glm-models-or-the-residual-d

To avoid any misunderstanding, can you explain how is it that two users have the same dataset and exactly the same code? Is it a school assignment?

Hi @StupidWolf They are some results from a group project yes, I guess someone else in the group is asking the same question, I will find out who and communicate with them to keep you from answering both of us!

@StupidWolf I see you say it indicates some of the terms are linear correlated, how come the glm function is not throwing up an error then?

Not linearly correlated, but linear combinations. It means that you don't have the samples for you to estimate certain effects. For example if all samples from origin A are also of substrate 1, then you cannot distinguish A and 1 right? https://stats.stackexchange.com/questions/35071/what-is-rank-deficiency-and-how-to-deal-with-it

@StupidWolf I think I get what you mean but that is as are as I am aware, not the case in this data. There might however examples of substrates only found in origin A or B, Could this cause the same effect?

@StupidWolf I just checked and there are 9 substrates found in origin B that are not in origin B

@StupidWolf I am artificially forcing this though but just putting a count of each substrate for each sample. and the samples are from either origin A or B

@StupidWolf Is there some way I can share the data with you?
If the substrate did not exist in that sample it was simply given a count of 0, so i did fill for missing counts (I think?) if I am understanding correctly.

Let us continue this discussion in chat.

When you fit a model, you need to check whether it's overly determined, meaning some terms cannot be determined. We can look at one of your samples, for example only for RB5 and Nocardia_nova_SH22a_RB20_RB56:
df = structure(list(X = c(226L, 227L, 228L, 229L, 230L, 231L, 232L,
233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 242L, 243L,
244L, 245L, 246L, 247L, 248L, 249L, 250L, 251L, 252L, 253L, 254L,
255L, 256L, 257L, 258L, 259L, 260L, 261L, 262L, 263L, 264L, 265L,
266L, 267L, 268L, 269L, 270L, 1126L, 1127L, 1128L, 1129L, 1130L,
1131L, 1132L, 1133L, 1134L, 1135L, 1136L, 1137L, 1138L, 1139L,
1140L, 1141L, 1142L, 1143L, 1144L, 1145L, 1146L, 1147L, 1148L,
1149L, 1150L, 1151L, 1152L, 1153L, 1154L, 1155L, 1156L, 1157L,
1158L, 1159L, 1160L, 1161L, 1162L, 1163L, 1164L, 1165L, 1166L,
1167L, 1168L, 1169L, 1170L), sample = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Nocardia_nova_SH22a_RB20_RB56",
"RB56"), class = "factor"), substrate = structure(c(10L, 37L,
29L, 11L, 24L, 9L, 39L, 30L, 42L, 14L, 21L, 40L, 38L, 45L, 28L,
5L, 22L, 3L, 41L, 16L, 32L, 27L, 18L, 26L, 20L, 34L, 7L, 12L,
23L, 17L, 44L, 35L, 4L, 19L, 25L, 13L, 1L, 33L, 31L, 43L, 36L,
6L, 2L, 8L, 15L, 10L, 37L, 29L, 11L, 24L, 9L, 39L, 30L, 42L,
14L, 21L, 40L, 38L, 45L, 28L, 5L, 22L, 3L, 41L, 16L, 32L, 27L,
18L, 26L, 20L, 34L, 7L, 12L, 23L, 17L, 44L, 35L, 4L, 19L, 25L,
13L, 1L, 33L, 31L, 43L, 36L, 6L, 2L, 8L, 15L), .Label = c("acylsphingosine",
"agarose", "alcohol", "alginate", "alpha-glucan", "amylaceous polysaccharides",
"arabinogalactans", "aribinose", "beta-glucan", "cellulose",
"chitin", "cutin", "cyclomaltodextrin", "D-galactose", "D-glucan",
"D-glucose", "dermatan sulfate", "dextran", "fructans", "fructose",
"galactan", "galactose", "galacturonide", "glucan", "glucose",
"glucuronoside", "glycan", "glycans", "glycogen", "hemicellulose",
"hyaluronate", "L-rhamnose", "lacto-N-tetraose", "levan", "lignin",
"maltose", "mannan", "mannose", "pectin", "peptidoglycan", "pullulan",
"trehalose", "Vanillyl-alcohol", "xanthan", "xylose"), class = "factor"),
    count = c(2L, 0L, 2L, 2L, 0L, 1L, 2L, 2L, 3L, 0L, 0L, 1L,
    0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 3L, 1L, 3L, 3L, 0L, 3L, 3L, 5L, 2L, 0L, 0L, 4L,
    0L, 0L, 1L, 2L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
    0L, 0L, 0L), origin = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("free", "isolate"
    ), class = "factor")), row.names = c(226L, 227L, 228L, 229L,
230L, 231L, 232L, 233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L,
241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L, 250L, 251L,
252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 260L, 261L, 262L,
263L, 264L, 265L, 266L, 267L, 268L, 269L, 270L, 1126L, 1127L,
1128L, 1129L, 1130L, 1131L, 1132L, 1133L, 1134L, 1135L, 1136L,
1137L, 1138L, 1139L, 1140L, 1141L, 1142L, 1143L, 1144L, 1145L,
1146L, 1147L, 1148L, 1149L, 1150L, 1151L, 1152L, 1153L, 1154L,
1155L, 1156L, 1157L, 1158L, 1159L, 1160L, 1161L, 1162L, 1163L,
1164L, 1165L, 1166L, 1167L, 1168L, 1169L, 1170L), class = "data.frame")

Now we look at whether how it is distributed:
df$sample = droplevels(df$sample)
table(df$origin,df$sample)

          Nocardia_nova_SH22a_RB20_RB56 RB56
  free                               45    0
  isolate                             0   45

You can see, the samples are either free or isolate, and never both. So you cannot discern the effect of origin from samples. This means your model is over-determined and hence you see the error above.
To know the effect of origin, you have to fit, i used a poisson for this example because if i use only the subset above, the counts are not overdispersed hence glm.nb throws some error:
fit = glm(count~ origin+substrate,data=df,family=poisson)
summary(fit)

 head(coefficients(summary(fit)))
                           Estimate   Std. Error       z value   Pr(>|z|)
(Intercept)           -2.266968e+01 2.930556e+04 -7.735624e-04 0.99938279
originisolate          6.931472e-01 2.970443e-01  2.333481e+00 0.01962291
substrateagarose       2.282092e-06 4.144430e+04  5.506408e-11 1.00000000
substratealcohol       2.280854e-06 4.144430e+04  5.503420e-11 1.00000000
substratealginate      2.278865e-06 4.144430e+04  5.498621e-11 1.00000000
substratealpha-glucan  2.295736e+01 2.930556e+04  7.833790e-04 0.99937495

You can see that in this model, free is the reference, and isolate is compared to free in the "originisolate". The coefficient 6.931472e-01 means isolate samples are higher than free samples, by exp(6.931472e-01) times.


text:
Bacterial DNA at the tail of transcriptome reads. What does that mean?
I am assembling a transcriptome obtained from the Internet. The transcriptome was extracted from a human cancer tissue that had been previously grafted into a mouse. I have detected that many assembled RNA sequences do not match precisely to their reference but have additional sequences at one of their tails containing bacterial DNA. It's not that there are reads pertaining to that bacterium (always the same one) but that the RNA sequences of the cell have incorporated fragments of bacterial sequences (the reads corroborate these insertions with high coverage), therefore I can't just filter out bacterial reads as there are not reads exclusively pertaining to bacteria but hybrids.
My question is how can it be possible? What possible mechanisms could have brought these bacterial fragments to the RNA of the cell. I have read in some papers that there is speculation that bacteria could be involved in the genesis of cancer cells. That would be an explanation if it wasn't because there are also mouse RNA for which the same circumstance applies: they also have bacterial sequences at the tails. I have also read in papers that there is speculation with the possibility that cancer cells are more susceptible to get DNA from bacteria transferred to their DNA. I am speculating that is the reason in this case (and somehow that weakness was transferred to contiguous mouse cells) but I would like to know if there is another possibility I am missing.
Thanks.

Sorry, I mean there are bacterial sequences at the tail of the contigs, not only the reads.

Sounds like an opportunistic infection, which bacteria? Just blast it, anything pathogenic, or just coliforms ... if its mixed  bits of genome from the same species - thats your answer. The mouse is very sick these things happen

The bacteria is Ralstonia solanacearum. It's mixed with mouse RNA and with human cancer RNA. there is no mix between mouse and human RNA. I'm trying to understand the mechanism by which the bacterial sequence ended up in the cancer and mouse RNA.

Usually it is like @swbarnes2 mentioned, adapter sequence. Can you check the fragment length of your aligned sequence? If they are short, then i think you can do an adaptor trimming first

@StupidWolf, the length of the reads is 140 bps and the contig is around 1400 bps for a transcript that should be 1280 bps long (GAPDH). That's only one of the transcripts. Not sure if you asked that. If I remove the adapters there are still more than 100 bps of bacterial sequence.

I was referring to the fragment length, not the just the read length. So you mentioned you downloaded it from ????internet????, so it helps to provide information on that. Is this a paired end library? If it is, you can get the fragment length from it and you will know whether you are sequencing into the adaptors

From time to time these situations are encountered, notably when PCR is being deployed in rare DNA samples.
In this instance the wiki states,

Ralstonia solanacearum is an aerobic non-spore-forming,
Gram-negative, plant pathogenic bacterium.".

I'd leave it. The simple explanation is there has been a bad mix up. The best explanation is the people loading the Illumina machine, designating the adaptors etc ... In this scenario it has nothing to do with bacterial genome integration, plant bacteria don't touch mice and aerobic bacteria wouldn't grow in cancer cells.
The second issue is what is a soil bacteria doing in a mouse cancer biological facility? For the staff this wouldn't be good news.
The singular biological  possibility is its a Pseudomonas  and some species are serious opportunistic pathogens of immune-weakened humans, but there are a lot of Pseudomonas species and this one doesn't ring any bells, BUT its mice so who knows? One species in particular is associated with hospital outbreaks, but a biomedical unit that isn't quite the same thing.

The pathogen in humans is called Pseudomonas aeruginosa. If hospitals get this infection it is BAD NEWS, for example burns victims are really susceptible and without intervention it kills them. If you have evidence that the bacterial equivalent in mice is the bacterial pathogen you have found that is a substantial step. Biological Safety Facilities (BSFs) DEFINITELY get infected. They hate it when this happens and can be very picky for fear of a contamination, but it certainly happens. If you have evidence of a major contaminent pathogen in the long run it would be of HUGE help to BSFs globally. No-one really bothers with the disease pathology of mice - because their mice and spend their entire lives going 'squeek' - but when they are models for cancer their status changes beyond recognition. So I'd erase my initial thoughts and suggest, what you have found is publishable and the application is BSF contamination - don't say that though say this could be the mouse equivalent of P. aeruginosa. Its certainly a publication and goodluck.

Thanks Michael. I guess your answer is correct but I am still trying to understand how a bad mix up can result in that data.

Hi @juanjo75es the best way to thank me is to upvote or even accept :-)  There is a possibility that the BSF (Biological Safety Facility) has been contaminated, it certainly happens and what you are looking at is the mouse equivalent of Pseudomonas aeruginosa (hospitals hate this disease).. If you can suggest this in publication that is the mouse equivalent of P. aeruginosa in humans for mice that is a major step forward. The BSFs of this world will not want to hear this news ... but who knows right? I'll update the post to incorporate this info

Thanks again, Michael! I will accept an answer once I really understand what's going on. Sorry but I am just developing software and have little experience in biology so I may be a bit slow.
I see now @swbarnes2 answer saying that this is not bacterial contamination (if I understand well) so ... I am still confused.

Thanks @juanjo75es

Looks like the sequence is just an adapter that was incorporated into the Ralstonia assembly  (and other ones) and that's why it matched bacterial genome. I guess I will have to accept another answer despite your kind help which is mysteriously uncommon in bioinformatics forums. I apreciate your interest to help.

Sure @juanjo75es ... if all the sequence you found was identical then the adapter is suspected. It was interesting to learn that they used a Pseudomonas type adapter from a rare plant pathogen.

What exactly is the sequence?  Are you sure it's not an artifact of the library prep?
5 seconds of googling turns up this:
https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/experiment-design/illumina-adapter-sequences-1000000002694-13.pdf
This sequence:
ACCACCGAGATCTACA
is on the adapters.
So is this
GACGGCATACGAGAT
Bacterial contamination would not put bacterial sequence at the end of host transcripts!  But adapter contamination sure could be.
Here's the index 2 adapter
AATGATACGGCGACCACCGAGATCTACAC[i5]TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Here's your sequence
CACCACCGAGATCTACACCTAAGCCTTCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTATCAACGCAGAGTACGGG
Here's it's blast against nr
ACCACCGAGATCTACACCTAAGCCTTCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTATC
It's all adapter.
Whoever did the transcriptome assembly didn't trim for adapters first.  That's all.

I am not sure of anything. I just want to know what are the possibilities. The data was downloaded from the Internet.
It's not always the same sequence. For example, I have this two at an extreme of a GAPDH RNA: CACCACCGAGATCTACACCTAAGCCTTCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTATCAACGCAGAGTACGGG and GACGGCATACGAGATCAGCCTCGGTCTCGTGGGCTCGGAGATTTTGT

Apart from the adapters, how do you explain the presence of the rest of the bacterial sequence?

Sorry swbarnes2 but I do not follow you. If I make BLAST search of that sequence it points to Ralstonia solanacearum with e-value 6e-23, query cover 83% and percentage of identity 100%. I made the BLAST search at  https://blast.ncbi.nlm.nih.gov/Blast.cgi
I certainly do not trim anything as I do not need that for assembling the transcriptome but that seems a more than reasonable answer. I just guess now why such long adapter and why using an adapter that matches bacterial dna. Hos do you make that blast search?

Could it be that the assemblies that BLAST finds wrongly include adapters in the assembled sequences? The sequence is repeated a lot of times in the assembly and in other organisms.

You should trim adapters away.  Looks like the Ralstonia sequencing group did not.

Could you please modify your answer pointing out that the problem is that the Ralstonia sequencing group did misassembly it (because of not trimming or for whatever other reason)? I'd like to accept your answer as it helped to find the real issue but the trimming thing seems just  circumstancial.


text:
Chromatin accessibility level on promoters with different CpG ratio calculation
I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).
So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.
In term of tool i have used chipseeker to annotate and do the above step.
So now if i have to find the CpG density as i have mentioned above how do i get that?

you can either calculate the GC content? And you can also check overlap with annotated cpg island from UCSC... Not very sure what is the question

can i get the GC content from  atac seq data ? if yes..how that is my question I mean what input do i require

did i get you wrong? It seems like you have annotated peaks. It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs (sorry I typed too quickly previously, it should be CGs instead of GC content)

"It seems like you have annotated peaks." yes correct. "It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs" okay that looks straight forward though can it be done in R ? can you suggest me some library

Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:
library(rtracklayer)
library(BSgenome.Hsapiens.UCSC.hg38)

peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"

extraCols <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1   10049-10326      * |
  [2]     chr1 180686-181040      * |
  [3]     chr1 186648-186903      * |
  [4]     chr1 191293-191604      * |
  [5]     chr1 267872-268146      * |
  [6]     chr1 586110-586324      * |
                                                   name     score signalValue
                                            <character> <numeric>   <numeric>
  [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
  [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
  [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
  [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
  [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
  [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035

Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:
library(Biostrings)
SEQ = getSeq(Hsapiens,gr)
head(SEQ)

A DNAStringSet instance of length 6
    width seq
[1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
[2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
[3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
[4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
[5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
[6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA

This is the CG, and you can decide how you want to bin your regions:
freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
hist(freq[,"CG"],br=100)

thank you veryyyyy much..i did look into couple of papers but it didn't occur to me how they did it..i will try this and update you..

"This is the CG, and you can decide how you want to bin your regions:" very naive  question the sequence also contain the "AT " how do you define CG?

sorry I typed too quickly, so in the code above, the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG, which is what you need, if you are interested in CpG?

yes i need CG for example this "TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC"  "the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG" this sequence which i have put also contains AT my question is is the function "dinucleotideFrequency" only look for CG only??

of course not.. see when i wrote this freq = dinucleotideFrequency(SEQ,as.prob=TRUE), freq is a matrix. the next line i called out freq[,"CG"], did you look at the matrix?

" freq[,"CG"]" im extremely sorry not sure how did i miss that!!!!  now i got it and big thank to make it clear to me...


text:
What is the best way to get a large number of RNA seq data from SRA in Python without being denied access
I have the following code to download data from the SRA using multithreading in Python. After running this a few times now (for testing purposes), I keep getting denied access to the data. Not sure how to fix this. In particular, cating the output files gives:
<?xml version="1.0" encoding="UTF-8"?><Error>AccessDenied<Message>Access Denied</Message><RequestId>cb0f0e98-cafb-1dd7-9b7b-d8c49756ec52</RequestId><HostId>QeDGVwBXYp61J0B4_OUTn7UsEsiQEec0n18DAeR0kaE</HostId></Error>

Here is my code:
import threading
import requests
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, url):
        threading.Thread.__init__(self)
        self.url = url
        self.result = None
        self.filename = url.split('/')[-1]
    def run(self):
        res = requests.get(self.url)
        with open(self.filename, 'wb') as f:
            f.write(res.content)
urls = [
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.2',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.2']

threads = [MyThread(url, ) for url in urls]
for thread in threads:
    thread.start()
for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

I once needed to download more than 1000 files from SRA, and I found a blog post that is very useful (https://reneshbedre.github.io/blog/fqutil.html). The idea is to use aspera to fetch SRA files. This is fast becasue of aspera. Then use fasterq-dump to convert SRA file to fastq. fasterq-dump is faster than fastq-dump partly because it does not compress the output file. I think compressing is the rate-limiting step here. After this, you can use pigz, which is a multi-threading version of gzip to compress the fastq files. Note that if you want to download files to a computing cluster, you can submit jobs for fasterq-dump and pigz. This makes things much faster. You cannot submit jobs for fetching SRA, because you need to connect to internet.

Thanks for sharing this. So I tried running prefetch with aspera as explained in the blog, but I get unrecognized option when I include the '-a' option. I do have aspera installed. Did you have this issue by any chance?

seems like aspera is no longer supported https://github.com/ncbi/sra-tools/issues/255

however, prefetch is still very fast, thanks for sharing this, update your answer to reflect this and I'll accept it!


text:
Trying to download a large number of files from ENA programmatically
I am downloading a large number of files from the ENA in Python using Multithreading. As an example, to be more concrete, I tried downloading 12 files simultaneously using 12 threads. I noticed that 10 files download relatively quickly, while 2 lag behind. Even when the 10 files are done downloading, the speed at which the other 2 files download is very slow. If I then try downloading these 2 files again, the speed will be a lot better. So what exactly is happening and how can I fix this?
PS: I tried decreasing the thread count to 8, so 8 files were downloading simultaneously, and I had 1 file that was lagging the rest again.
PPS: I then tried decreasing to 6 threads, and again the same issue. One file will lag. I don't have this issue if I directly download from Firefox, and I get a significant speed up, but I'm trying to automate the process so downloading from Firefox isn't an option.
For reference, here's the code:
import queue
import threading
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, name):
        threading.Thread.__init__(self)
        self.name = name
    def run(self):
        print('Starting thread %s.' % self.name)
        process_queue()
        print('Exiting thread %s.' % self.name)

def process_queue():
    while True:
        try:
            url = my_queue.get(block=False)
            filename = url.split('/')[-1]
            with closing(request.urlopen(url)) as r:
                with open(filename, 'wb') as f:
                    shutil.copyfileobj(r, f)
        except queue.Empty:
            return

# setting up variables
urls = [
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/001/SRR8240861/SRR8240861_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/001/SRR8240861/SRR8240861_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/002/SRR8240862/SRR8240862_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/002/SRR8240862/SRR8240862_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/003/SRR8240863/SRR8240863_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/003/SRR8240863/SRR8240863_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/004/SRR8240864/SRR8240864_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/004/SRR8240864/SRR8240864_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/005/SRR8240865/SRR8240865_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/005/SRR8240865/SRR8240865_2.fastq.gz',
        ]

# filling the queue
my_queue = queue.Queue()

for url in urls:
    my_queue.put(url)

# initializing and starting num_threads threads
num_threads = 8
threads = []

for i in range(num_threads):
    thread = MyThread(i)
    threads.append(thread)

for thread in threads:
    thread.start()

for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

Their web server may be throttling simultaneous connections. Generally, multiple incoming requests from one client is not good.

@burger I think you're right. I added a time.sleep(1) after thread.start() and it seems like it's better now. Though I'll have to test some more

@burger hmm so I performed some tests and it still seems like I have an issue. This is despite the fact that I have time.sleep(60) Firefox is still a lot faster at downloading the files than ftp


text:
How can I obtain FTP links to studies in ENA given a BioProject ID?
How can I programmatically obtain ftp links to RNA seq fastq files in ENA? Here's an example of a link that I would be interested in obtaining:
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz
In particular, is there some tool that, given the BioProject ID (here, PRJNA506829), would be able give me all ftp links for the runs in the project, or would I need to write a web scraper to do it?

Here are some command lines I used for that purpose in bash.
Simply prepare a text file containing each accession number (SRR/ERR) you want and create a for loop. Here I used prozilla to speed up downloads but you may use wget either.
for index in $(cat list_of_accessions) ; do

if [ ${#index} -eq 9 ]; then
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_1.fastq.gz
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_2.fastq.gz
            else
            if [ ${#index} -eq 10 ]; then
                    proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_1.fastq.gz
                            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_2.fastq.gz
                    else
                    if [ ${#index} -eq 11 ]; then
                     proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_1.fastq.gz
                                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_2.fastq.gz
                    else
                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_1.fastq.gz
                                                proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_2.fastq.gz
                        fi
                fi

            fi
done
```

They have an API you can interact with.
If you need to get files for only a few different projects:
Search for your accession ID in the browser (leading to https://www.ebi.ac.uk/ena/browser/view/PRJNA506829). Filter the show selected columns to only fastq_ftp, click download tsv to get the list of ftp links.

pysradb can fetch ENA/SRA fastq/bam links (if available):
$ pysradb metadata SRR8240860 --detailed
run_accession study_accession experiment_accession experiment_title                                               experiment_desc                                                organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument   total_spots total_size  run_total_spots run_total_bases run_alias      sra_url_alt1                                    sra_url_alt2                                    sra_url                                                                                experiment_alias source_name strain/genotype developmental stage ena_fastq_http ena_fastq_http_1                                                                 ena_fastq_http_2                                                                 ena_fastq_ftp ena_fastq_ftp_1                                                                     ena_fastq_ftp_2
SRR8240860    SRP170618       SRX5059122           GSM3487689: fer-15 Day1 rep1; Caenorhabditis elegans; RNA-Seq  GSM3487689: fer-15 Day1 rep1; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              PAIRED         SRS4075216       N/A          HiSeq X Ten  40494415    4755049874  40494415        12148324500     GSM3487689_r1  gs://sra-pub-src-3/SRR8240860/RRA0719_R2.fq.gz  s3://sra-pub-src-3/SRR8240860/RRA0719_R2.fq.gz  https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR8240860/SRR8240860.1  GSM3487689       whole worm  fer-15(b26ts)   Adult day 1         N/A            http://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz  http://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz  N/A           era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz  era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz

You can do this in one line using NCBI's EntrezDirect e-utils:
esearch -db sra -query PRJNA506829 | efetch -format runinfo
The output will be in csv, but you can format it to tab using csvkit's csvformat -T


text:
Calculating Isoelectric point from a multifasta file
I need to calculate Isoelectric point from a multifasta file, is there any python code or web tool that allows me to do that?

Can you use Python (biopython has a very straightforward method for pI calculations, Bio.SeqUtils.IsoelectricPoint)?

You can use EMBOSS Pepstats for this. It takes a multi-fasta as input and produces a table that contains various statistics related to the protein, including isoelectric point.


text:
How to fragment genomes into non-overlapping sequences of differing sizes?
I have bacterial chromosomes/plasmids that I have downloaded off refseq in a multi-fasta that I would like to fragment into non-overlapping fragments of between 1 kb and 15 kb.
I have tried googling for related tools but to no avail. Does anyone know how I would achieve this? Something very similar is done in this paper:
https://link.springer.com/article/10.1007/s40484-019-0187-4
Thanks in advance!

You could use numpy.random.randint to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.
Here's some code:
#!/usr/bin/env python

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next < len(sequence):
        new_next = np.random.randint(low, high) + next
        if new_next > len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = ""
sequence = ""
for line in sys.stdin:
    if line.startswith('>'):
        if len(sequence) > 0:
            process_record(header, sequence, low, high)
            sequence = ""
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) > 0:
    process_record(header, sequence, low, high)

To use it:
$ ./13879.py 1000 15000 < in.fa > out.fa

Fragments will be disjoint (non-overlapping, contiguous) subsequences.
I add a check to determine if the last fragment of a FASTA record would be shorter than low (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:
if new_next > len(sequence) - high:
    new_next = len(sequence)

You could change this logic to discard the last fragment (i.e. break), or simply allow fragments less than 1kb in size, depending on what you need.
EDIT
Here's a slight modification that should allow to you preserve fragment sizes within your bounds.
A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into n fragments:
   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|

We may get unlucky and sample a random integer for fragment n-1 (frag(n-1)) which makes the length of frag(n) less than low (1 kb in your example).
In my first suggested solution, frag(n-1) and frag(n) are appended to one another. This is why you get some fragments longer than high (15 kb).
Another option is to discard frag(n), which is of length less than low (1 kb). Here's some code:
#!/usr/bin/env python

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

"""
Globals
"""
fragment_lengths = []

"""
Process a record
"""
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = ""
    while last < len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)
        last = next
        fragment += 1
    """
    Uncomment to compare input sequence with fragments
    """
    #sys.stderr.write('In > {}\nOut > {}\n'.format(sequence, fragment_sequence))

"""
Parse records
"""
def parse_records(in_fn, low, high):
    header = ""
    sequence = ""
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('>'):
                if len(sequence) > 0:
                    process_record(header, sequence, low, high)
                    sequence = ""
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) > 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    """
    Comment below if testing the distribution of fragments
    """
    parse_records(in_fn, low, high)
    """
    Uncomment to test the distribution of fragments from
    running fns on same input 100k times
    """
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == "__main__":
    main(in_fn, low, high)

If frag(n) is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then rejection sampling might be an option.
With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.
If you need that, follow up and I'll take another look.

Thanks for your answer! Do I need to remove all the fasta headers? I tried it on a multifasta and ended up with a few sequences over 15 kb. The distribution looks like this: https://imgur.com/a/FFr8xNv

No, you do not need to remove headers. You could end up with a few sequences over 15kb, if the last fragment cannot be sized between low and high (i.e. between 1 and 15kb). I'll propose a solution in the edit.

Thanks that worked perfectly!

Hi @StanleyHo please do sign the answer as "asscepted". I'll be making use of this code too, I like the args and kwargs

I've previously written my own script to fragment DNA sequences (input as either fastq or fasta files) into identical-length fragments, with a given overlap. I used this because there was a particular assembler that demanded equal-length sequences, and I was working with nanopore reads, but it's since been useful for other purposes as well.
Because it is most concerned about making lengths identical, there may be some forced overlap when setting the overlap to zero. The overlap will only be exactly zero when the sequence length is an exact multiple of the fragment size. Here's an example use:
zero overlap
normalise_seqlengths.pl -fraglength 25 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 6 fragments, 0 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGT
+
FFFFFFFFFFFFFFFFF:,FFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:25-50  1:N:0:CGCAAGAA+ACTCCTAA
GGGAGCTGTCTGGGGTAACCTGGCT
+
:FFFF:F:FFFFF::FFFFF,::FF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:50-75  1:N:0:CGCAAGAA+ACTCCTAA
CTAGGACCAGGTGGAAAAAGCACCT
+
F:FFFFFFFFFFFFFFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#3:75-100  1:N:0:CGCAAGAA+ACTCCTAA
TGGCTGGCCCAGCCCTGTTCTTCTC
+
:FF,FFFFF,:,FFF::F:FFF,:F
@A00187:180:HLGNGDSXX:3:1101:8486:1031#4:100-125  1:N:0:CGCAAGAA+ACTCCTAA
CTGGGGAGGGAAGGGCAGGTCCCCT
+
F:FFFFFFFFF,FFFFF:FFF,F,,
@A00187:180:HLGNGDSXX:3:1101:8486:1031#5:125-150  1:N:0:CGCAAGAA+ACTCCTAA
AATGCCTGAACCACATGGTGGGCCA
+
FFF,F:FFFFFFF:F,,,FFF,F,,

partial overlap
normalise_seqlengths.pl -fraglength 63 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 3 fragments, 20 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGTGGGAGCTGTCTGGGGTAACCTGGCTCTAGGACCAGGTG
+
FFFFFFFFFFFFFFFFF:,FFFFFF:FFFF:F:FFFFF::FFFFF,::FFF:FFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:43-106  1:N:0:CGCAAGAA+ACTCCTAA
CCTGGCTCTAGGACCAGGTGGAAAAAGCACCTTGGCTGGCCCAGCCCTGTTCTTCTCCTGGGG
+
FF,::FFF:FFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFF,:,FFF::F:FFF,:FF:FFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:87-150  1:N:0:CGCAAGAA+ACTCCTAA
CCCTGTTCTTCTCCTGGGGAGGGAAGGGCAGGTCCCCTAATGCCTGAACCACATGGTGGGCCA
+
FFF::F:FFF,:FF:FFFFFFFFF,FFFFF:FFF,F,,FFF,F:FFFFFFF:F,,,FFF,F,,


text:
Kraken2 > OTU format > Phyloseq
A collaborator has passed me over Kraken2 outputs *.report and *.kraken, from a metatranscriptomic sequencing experiment conducted on the minION.
I would like to make a tree if the data using a standard phylogenetics package such as phyloseq, however, I am having trouble making an alignment via the OTU table..
There has been a similar question on here already, but the author has resolved it and the solution is not available.
Does anyone have any pointers?

An OTU table is here https://github.com/smdabdoub/kraken-biom, what you require is an alignment for phyloseq


text:
Difference between "trans-ethnic" and "cross-ancestry"
A quick terminology question today. I see the terms "cross-ancestry" and "trans-ethnic" used seemingly interchangeably in literature. Is there any real difference between those two or are the 100% the same?


text:
downloaded software in Conda environment does not show in ls command
I installed gzip package in Conda environment, then  to make sure the software was there, I used "ls" command, it was not in the list. Tried to install again, system says all the packages have already installed, Why? How do I check the list of software installed in the environment? Thanks!

Just type
conda list gzip

OR
conda list | grep gzip

OR if its installed in a given environment
conda activate myenv

or else just activate myenv system depending followed by either of the first two commands.
It is better to use conda-navigate a very nice and under used GUI, and would work in this particular instance. The details are found in an earlier post I made on esentially the same question here together with a pretty picture.
ls will only work if you type ls ~/Anaconda3/bin | grep gzip , but depends on the installation, e.g. you might have installed miniconda.
Having said all that conda is very cool, but whatever you do NEVER use pip or pip3 from this moment onwards EVER, it messes up your conda and you'll need to env your way out of it.


text:
nanopore QC measures on fastq file
I am working with virus nanopore sequencing data. I have to remove host genome from the reads. I want to run the QC on the fastq file after host genome subtraction. I have used nanoplot. It is slow. I dont have summary.txt. Is there any nanopore read specific qc package which can directly take nanopore fastq files.

Oh yeh, NanoPlot can be pretty slow when you only have fastq data. I'm happy to look at contributions on how to make it faster. Do you have a machine with multiple cores? Because if your data is split across multiple fastq files then NanoPlot can read those in parallel.

I have 6 core with 12 threads. How can I split the fastq files and run them in parallel?

Splitting can be done e.g. with https://fastqsplitter.readthedocs.io/en/stable/, or the Unix split tool.


text:
How to use pyopenms.FeatureFindingMetabo() for peak picking in pyOpenMS?
I found this documentation on pyOpenMS mass spec library. I want to try different peak picking algorithms. Unfortunately, the docs are not very detailed and don't show how to use the other peak picking algorithms and the codebase itself is also quite cryptic for me.
Any idea how to use, for example, the FeatureFindingMetabo algorithm insead of the centroided one?
from urllib.request import urlretrieve
# from urllib import urlretrieve  # use this code for Python 2.x
gh = "https://raw.githubusercontent.com/OpenMS/OpenMS/develop"
urlretrieve (gh +"/src/tests/topp/FeatureFinderCentroided_1_input.mzML", "feature_test.mzML")

from pyopenms import *

# Prepare data loading (save memory by only
# loading MS1 spectra into memory)
options = PeakFileOptions()
options.setMSLevels([1])
fh = MzMLFile()
fh.setOptions(options)

# Load data
input_map = MSExperiment()
fh.load("feature_test.mzML", input_map)
input_map.updateRanges()

ff = FeatureFinder()
ff.setLogType(LogType.CMD)

# Run the feature finder
name = "centroided"
features = FeatureMap()
seeds = FeatureMap()
params = FeatureFinder().getParameters(name)
ff.run(name, input_map, features, params, seeds)

features.setUniqueIds()
fh = FeatureXMLFile()
fh.store("output.featureXML", features)
print("Found", features.size(), "features")

Hi Sören, We build a Python pipeline with FeatureFinderMetabo, it was like a year ago so I can not recall the details right now, but you can see the entire thing here: https://github.com/saezlab/lipyd/blob/master/src/lipyd/msproc.py
If you will be still having issues, maybe in the weekend I can check it again.

This is a draft answer, if the OP still need needs I will extend it later.
Me too I found the documentation and examples of pyopenms very scarce. At the same time  OpenMS is huge, sometimes you find key information not in various other sources instead of the official documentation. This is especially true if it comes to metabolomics.
We built a Python OpenMS lipidomics LC MS/MS preprocessing pipeline, contained in this file:
https://github.com/saezlab/lipyd/blob/master/src/lipyd/msproc.py
The top class is the one called MSPreprocess, it implements the workflow, the other classes implement its steps.
Briefly, you start from profile mode mzML files, then

Peak picking with PeakPickerHiRes results centroided data
FeatureFindingMetabo assembles the traces in a metabolomics
compatible way
A MapAlignment algorithm alignes the traces across
scans
A FeatureGrouping algorithm alignes the features across multiple experiments

Our software mentioned above (lipyd) outputs data frames from the m/z, RT and intensity.
About OpenMS I can recommend to read the C++ code and docs. The Python classes and methods are automatically built minimal wrappers, everything works the same way as in C++.

Hi, thank you for sharing the link. I will have a look at it.  Do you mind sharing what version of pyopenms was used here?

I only know both openms and pyopenms on my computer are 2.5.0, and I think I haven't updated them since last summer

It's hard to understand what is going on and when FeatureFinderMetabo is called and with which arguments. As it is so encapsulated in the method_key wrapper.

Anyway, it was super helpful. I composed my own answer out of your code, which is more concise. You could comment on that in case I made a mistake or forgot something.

My minimal procedure based on the code that @deeenes shared for a single file. Comments to improve this solution are highly appreciated:
files = list_of_filenames

feature_map = oms.FeatureMap()
mass_traces = []
mass_traces_split = []
mass_traces_filtered  = []

exp = oms.MSExperiment()

options = oms.PeakFileOptions()
options.setMSLevels([1])

fh = oms.MzXMLFile()
fh.setOptions(options)

# Peak map
peak_map = oms.PeakMap()

file = files[0]

print('# Filename:', basename(file))
fh.load(file, exp)

print('# Spectra:', len( exp.getSpectra() ))
print('# Chromatograms:', len( exp.getChromatograms() ) )

for chrom in exp.getChromatograms():
    peak_map.addChrom(chrom)

for spec in exp.getSpectra():
    peak_map.addSpectrum(spec)

mass_trace_detect = oms.MassTraceDetection()
mass_trace_detect.run(peak_map, mass_traces, 100000)

print('# Mass traces:', len(mass_traces) )

elution_peak_detection = oms.ElutionPeakDetection()
elution_peak_detection.detectPeaks(mass_traces, mass_traces_split)
print('# Mass traces split:', len(mass_traces_split) )

feature_finding_metabo = oms.FeatureFindingMetabo()
feature_finding_metabo.run(
            mass_traces_split,
            feature_map,
            mass_traces_filtered)

print('# Mass traces filtered:', len(mass_traces_filtered) )
print('# Features:', feature_map.size() )

oms.FeatureXMLFile().store('FeatureFindingMetabo.featureXML', feature_map)

>>>
# Filename: CA_1.mzXML
# Spectra: 989
# Chromatograms: 0
# Mass traces: 18476
# Mass traces split: 25194
# Mass traces filtered: 0
# Features: 23549


text:
Removing data from OrthoFinder output prefix
I am using OrthoFinder to find orthogroups between species. I am using the latest version unpacked through tar as per their instructions.
Orthofinder creates an output folder according to the current date, ie, Results_Dec25. However,
I want to integrate orthofinder in snakemake, so having a stable output directory is necessary.
The -o option allows specification of an output prefix, but the dated folder is then created within.
Is it possible to force a given name to the output folder?
Or is there a way in snakemake to allow ambiguous output directory names?
I know this would be possible via editing the source code of OrthoFinder, but this workflow needs to be reproducible without faff.
Many thanks
This question has already been asked here but remains unanswered.
Is there an option in Orthofinder to choose the output directory name?

The latest version solves the issue.
Source : https://github.com/davidemms/OrthoFinder/issues/403
Latest Version: https://github.com/davidemms/OrthoFinder/releases/tag/2.4.0


text:
How to install/import Biopython into Python 3.8/ PyCharm IDE
I'm new to Python/Biopython programming. I'd like to install Biopython into my PyCharm IDE but I've tried all instructions on this, but it's not confirming it. I have the available biopython file already extracted. There's a SyntaxError for 'install' in:
pip install biopython
What's the proper way of installing and importing biopython? And how does one confirm the installation.

You probably should mention your operating system.

Could you please show the complete error text? Which python version do you use?

A SyntaxError is an exception raised by python, so I guess you ran pip within the python interpreter rather than from your shell/terminal

@WouterDeCoster It's Windows 10 Pro, 64-bit OS, 4gb RAM

@Mr_Z Python 3.8, PyCharm 2020.1.2.        Here: pip install biopython
File "", line 1
pip install biopython
^
SyntaxError: invalid syntax

Install Anaconda3 the instructions are here. Then making sure you are connected to the internet type
conda install -c anaconda biopython

Once Anaconda is installed don't attempt to use pip, if you do then you need to set an 'environment' in Anaconda and work from that environment.
I suspect you used pip or pip3 via
pip install biopython


text:
maftools tcgaCompare plot: x-axis labels getting clipped
I have a {maftools} plot:
whatever I try to save that completely I see the down part of plot is trimmed like

I have tried both save as pdf and save as image options which did not differ
Any suggestion please?
Plot comes by this code
tcgaCompare(maf = laml,cohortName = "Cohort",rm_zero = FALSE,decreasing = TRUE,logscale = TRUE,primarySite = TRUE,col = c("grey30", "red"),medianCol = "black")

One way is to increase the outer margins of plot.
par(oma=c(7,0,0,0))
tcgaCompare(maf = laml, primarySite = TRUE)


text:
Kraken2 or metaphlan2 report to phyloseq
I wish to convert a kraken2 or metaphlan2 report into the R package to analyze using Phyloseq.
I have performed obtained an importation and Phyloseq analysis, but the result was not correct.
library("phyloseq")
library("biomformat")
otu_Rat11 <- read_biom("Rat11.biom")

otumat_Rat11 = as(biom_data(otu_Rat11), "matrix")
OTU_Rat11 <- otu_table(otumat_Rat11, taxa_are_rows=TRUE)
str(OTU_Rat11)
#Formal class 'otu_table' [package "phyloseq"] with 2
#slots ..@ .Data : num [1:1018, 1] 942 433 3007 116765 15103 ... .. ..-
#attr(*, "dimnames")=List of 2 .. .. ..:ℎ[1:1018]"186826""33958""1578"
#"33959".........:chr[1:1018] "186826" "33958" "1578" "33959"......... :
# chr "sa1" ..@ taxa_are_rows: logi TRUE

taxmat_Rat11 = as.matrix(observation_metadata(otu_Rat11), rownames.force=TRUE)
TAX_Rat11 = tax_table(taxmat_Rat11)

physeq_Rat11 = phyloseq(OTU_Rat11, TAX_Rat11)
str(physeq_Rat11)
#Formal class 'phyloseq' [package "phyloseq"] with 5 slots ..@ otu_table:Formal
#class 'otu_table' [package "phyloseq"] with 2 slots .. .. ..@ .Data : num
#[1:1018, 1] 942 433 3007 116765 15103 ... .. .. .. ..- attr(*, "dimnames")=List
# of 2 .. .. .. .. ..:ℎ[1:1018]"186826""33958""1578""33959".............: chr
#[1:7] "taxonomy1" "taxonomy2" "taxonomy3" "taxonomy4" ... ..@ sam_data : NULL ..
#@ phy_tree : NULL ..@ refseq : NULL –

Now I try to obtain @ sam_data
sampledata = sample_data(otu_Rat11)
#Error in access(object, "sam_data", errorIfNULL) : sam_data slot is empty.

What do with sample data and how to merge several phyloseq-format datasets from different metagenomes?

Is it possible to have an example file for Rat11.biom? It's easier to debug with a reproducible result.
What do you mean by "the result was not correct"?
What's your code for merging several phyloseq-format datasets from different metagenomes?

I'm working on this procedure with some difficult. Could you please tell me how do you convert kraken2 report files into the biom format?

I had received the same error. My metadata file was in a matrix format, I converted it into a data.frame and it worked. Pls try and see if it works.


text:
scRNA-Seq: Account for sequencing depth and gene length?
Task: Normalize a single-cell RNA-Seq dataset to account for sequencing depth and gene-length.
For UMI-count based protocols (like 10x) that don't suffer from gene-length biases, there are various methods developed to account for sequencing depth e.g. sctransform from Seurat.
For read-count based protocols (like Smart-Seq2), SCnorm can for example be used to normalize for sequencing depth. However, read-count based protocols can suffer from gene-length biases that are not accounted for in SCnorm. Furthermore, the method of SCnorm expects raw counts, so I don't know if I can length-normalize the counts prior to using the method. Technically, TPMs do account for both, but it is a global scale factor method, which has other drawbacks.
Question: How can I normalize my read-count based datasets to account for gene-length biases and sequencing depth without using a global scale factor method?


text:
Download multiple fastq files using fastq-dump
I want to download the following fastq files at the same time in Salmon:
 - SRR10611214
 - SRR10611215
 - SRR10611215
 - SRR10611216
 - SRR10611217

Is there a way to do this using a bash for loop or fastq-dump? or prefetch

https://www.biostars.org/p/347050/ you can refer to this answer. Here you are trying each read separately. it would be better if you can download and split the read.

A sample code is given in the salmon documentation as follows. Source
#!/bin/bash
mkdir data
cd data
for i in `seq 25 40`;
do
  mkdir DRR0161${i};
  cd DRR0161${i};
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz;
  cd ..;
done
cd ..

This could be modified as follows.
#!/bin/bash
mkdir data
cd data
for i in `seq 14 17`;
do
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_2.fastq.gz;
done
cd ..

You can save the code as a shell script and run it from the linux terminal. For example bash download_fastq.sh

Make a list.txt file containing a single column of SRA numbers to download.
then:
for i in $(cat list.txt); do echo $i; date; fasterq-dump -S $i; done
It works well to use NCBI's web interface to find SRA samples of interest, download and open findings in Excel, then copy single column containing SRA numbers and paste into list.txt using document editor such as vim.
After downloading including the "R" can be nice:
for i in *_1.fastq; do mv $i ${i%_1.fastq}_R1.fastq; done
for i in *_2.fastq; do mv $i ${i%_2.fastq}_R1.fastq; done
and zip:
pigz *fastq
If needed a conda option for downloading fasterq-dump:
conda install -c bioconda sra-tools

You can use parallel.
parallel -j 3 fastq-dump {} ::: SRR10611214 SRR10611215 SRR10611215 SRR10611216 SRR10611217

The option -j says how many jobs should maximal run parallel. So in this case maximal 3 identifier would be handled at the same time.
How many jobs you can run parallel depends on your machine.
You can also take a look at parallel-fastq-dump.


text:
P-value correction when evaluating correlation between gene and miRNA expression
First of all I apologize without the question is very basic, I am taking my first steps in bioinformatics.
Data information
We are evaluating the correlation (using the Pearson, Kendall or Spearman method) between gene expression and miRNA expression using the corAndPvalue function of WCGNA.
The resulting structure would be a DataFrame containing all combinations between each gene with each miRNAs, containing the following columns:
Gene     miRNA      Correlation P-value
Gen_1    miRNA_1    0,959       0.00311
Gen_1    miRNA_2    -0,039      0.1041
Gen_1    miRNA_3    -0,344      0.0021
Gen_2    miRNA_1    0,1333      0.00451
Gen_2    miRNA_2    0,877       0.07311
...

Question
Considering the huge number of correlation tests we are going to evaluate, we need to adjust the p-values to avoid correlations due to chance. Bonferroni does not seem to be the best solution, so we would use Benjamini-Hochberg method (BH). The question is:
The BH correction for the Gen_1 x miRNA_1 combination, should consider the p-values of all combinations that include Gen_1 (Option 1), or should consider all the p-values of all the genes x miRNA combinations (Option 2)?
For example, let's assume an expression dataset of 20,000 genes and another of 15,000 miRNAs
Option 1:
To adjust Gen_1 x miRNA_1 we would use 15,000 p-values (Gen_1 x miRNA_1, Gen_1 x miRNA_2, ..., Gen_1 x miRNA_15000).
Option 2:
To adjust Gen_1 x miRNA_1 we would use 300,000,000 p-values (Gen_1 x miRNA_1, Gen_1 x miRNA_2, ..., Gen_1 x miRNA_15000, Gen_2 x miRNA_1, Gen_2 x miRNA_2, ..., Gen_2 x miRNA_15000 and so on).
Suplementary question
Documentation of the method fdrcorrection from Python Statsmodels library suggests that for negative correlations (that could be frequent in a mRNA x miRNA correlation analysis) Benjamini-Yekutieli would work better; is that right? Or Benjamini-Hochberg method would be appropiated for this case?
Any kind of help would be much appreciated, thanks in advance!

There's a good answer in CrossValidated forum! Link

I made the same question in CrossValidated forum and got an excellent answer!
The important part:

You need to correct for all of the comparisons you are doing. So if
that's 300,000,000 comparisons you need to correct for that many
multiple comparisons.

For more information check the answer in the link above


text:
How do I get cromwell to create soft-links?
The file system I'm running on doesn't allow the creation of hard-links.  Cromwell is trying to create hard-links to directories, and when that fails, it's copying them.  Which is leading to disk quota errors.
Is there a way that I can simply tell cromwell "try soft links first"?  Without creating a whole config file?

Hi could you tell us what your disk quota is please?

@Michael 4.9 TB.  Which my pipeline fills up shortly before it's supposed to finish

It is possible to set configuration values via the command line. I haven't tried this, though (and it will probably make you command pretty long...)

@DavyCats we broke down and created a config file.  Problem is that for some reason cromwell insists on creating ".tmp" directories with complete copies of the working directory


text:
VCF spec: is it possible to have other alleles in addition to the MISSING value ('.') in ALTS?
I'm writing a VCF parser, so I have to consider and handle all corner cases regardless of how contrived they may seem. The specification is a bit unclear about the MISSING value ('.') in the ALTS column:

Options are base Strings made up of the bases A,C,G,T,N,*, (case insensitive) or a MISSING value ‘.’  (no variant) or an angle-bracketed ID String (“”) or a breakend replacement string.

I've seen examples with a single dot in the ALTS column:
4       31789170        PTV021  G       .       77      PASS    .

The question is whether the following data lines are also valid:
1       12345   ID1     A       .,T,.   22.88   PASS    .
1       12346   ID2     G       C,.     22.88   PASS    .

In other words, does MISSING indicate that the entire ALTS field is missing, or does it mean that there's a missing allele?
By extension, how do I represent the case when there's a single dot in the ALTS field (as in the first example)? Is it an empty list (because the whole field is MISSING) or is it a list containing a MISSING value? In other words, is it [] or ["."]?
At the moment, the parser does not validate the contents of the ALT field and does not even attempt to resolve any ALT IDs from the header. Interpretation of the ALT values is left entirely to the client code.

Is your parser intended to parse vcfs containing just SNPs, or different kinds of structural variants as well?

At the moment, the parser does not validate the contents of the ALT field and does not even attempt to resolve any ALT IDs from the header. Interpretation of the ALT values is left entirely to the client code.


text:
How can I restrict the search space when querying GEO?
I am using the following command to retrieve GEO Datasets entries for daf-2 c. elegans that have as their experiment type High Throughpout sequencing:
esearch -db gds -query "daf-2" | efilter -query "expression profiling by high throughput sequencing [DataSet Type] AND Caenorhabditis elegans [organism]"

However, among the results that I got, I noticed that I got this:
23. Impaired Insulin-/IGF1-Signaling Extends Life Span by Promoting Mitochondrial L-Proline Catabolism to Induce a Transient ROS-Signal
(Submitter supplied) Transcriptome profiling of three models with impaired insulin/IGF1 signaling. 1. Deep sequencing of endogenous mRNA from Caenorhabditis elegans N2 var. Bristol (wildtype) and daf-2(e1370) mutant; 2. Deep sequencing  of endogenous mRNA from murine embryonic fibroblasts (MEF)  wildtype and irs1-/- knockout; 3. Deep sequencing of endogenous mRNA from murine embryoinic fibroblast (MEF) insr+/- -lox and insr+/- knockout  Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de)
Organism:   Mus musculus; Caenorhabditis elegans
Type:       Expression profiling by high throughput sequencing
Platforms: GPL11002 GPL13776 14 Samples
FTP download: GEO (CSV) ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE36nnn/GSE36041/
SRA Run Selector: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA151765
Series      Accession: GSE36041 ID: 200036041

As you can see, in the Organism field, there are both Mus musculus and Caenorhabditis elegans. How can I restrict my search so that I only get entries that exclusively work on Caenorhabditis elegans? One naive way to do this would be to include NOT Mus Musculus in the query, but of course that would mean that I could get something else too.
Another way I can think of solving this is just to write a script that would do the extra filtering using regexes, but I was wondering if there's a simpler solution using the e-utilities functionalities

That particular study was done using both worm and mouse samples. If you are ultimately interested in SRA data, you should just pipe the output of this step to elink and get SRA run table. You can either add an organism filter after that or just subset the run table.

@vkkodali but why would filtering from SRA's side be different from filtering from GEO's side?

@vkkodali I noticed that it does filter out other organisms, but I just don't understand why filtering from SRA is different from GEO. Also, I'm filtering again by "Expression profiling by high throughput sequencing [DataSet Type]" because of the same issue as with organisms, but when I do this, I get 0 results. Any clue as to why this might be happening?

Filtering from SRA would filter it at the run level where as filtering in GEO at the series level will find any series that contain at least one sample that is from worm.

@vkkodali I see that makes sense. What about the other filter?

You can use the following query to first search GEO to get a list of GEO Series of interest, then find all linked SRA runs and perform a second filter step to keep only the SRA runs that satisfy another set of criteria as shown below.
esearch \
  -db gds \
  -query "daf-2[all fields] \
    AND expression profiling by high throughput sequencing [DataSet Type] \
    AND Caenorhabditis elegans [organism]" \
  | elink \
  -db gds -target sra \
  | esearch \
  -query "(#3) \
    AND Caenorhabditis elegans [organism]" \
  | efetch -format runinfo

Here, I am combining two independent queries using the (#3) term in the second query. More information about this is in the "Combining Independent Queries" subsection of "Searching and Filtering" section of the Entrez Direct documentation.

could you explain why you have a -db gds with elink?

It’s just a habit. If you are piping results from another edirect tool, you don’t need to specify the db. But if you are running elink as a stand-alone command then you need the db.


text:
Adding treatment groups via metadata to Seurat object?
I'd like to add metadata to 6 individual Seurat objects so that after I merge the objects into one, I can later label or split by using these identifiers. For example, I'd like to append an age group and then interval across these 6 objects. For adding an interval, I've tried using the below:

AddMetaData(control, metadata = 1hr, col.name = Interval)

However, this error will occur:
Error in .AddMetaData(object = object, metadata = metadata, col.name = col.name) :
  object '1hr' not found

Appreciate any help, thank you!

You should have fed "1h" to your function call but you fed 1h instead. The latter is looking for an object with that name and as the error shows, there is no such object.


text:
Human Herpes Virus, fusion motif in protein vs. SARS-CoV-2?
Is it possible that the fusion with cell propriety of Human Herpes virus, I see X14112.1 in GenBank is due to protein UL36 (tegument of virus) starting in linear position 2685 ???
H R R A R S L A

there are similarities with SarsCov2 virus protein S in name MN908947.3 in GenBank position 681.
P R R A R S V A

It is strange but the now (21-12-2020) the new Sars2 Great Britain strain is more near to that UL36 piece, it seems already:
H R R A R S V A

Herpes virus protein UL36 and SarsCov2 virus protein S[1200 position], have  too
L Q E L G K

i find in internet they call it "Glutaminase kidney isoform, mitochondrial"

If it's on a transmembrane protein and that is close to the transmembrane region that could be a di-arginine motif for ER/Golgi/cell surface translocation. It is super common.

You will always find similarities between sequences subject to evolution laws. New genes are not random sequences but sequences obtained reformatting other sequences, even from other species and with different functions. But these similarities have usually a very old origin. In the case of viruses, the fast mutation rates can make these sequences return to an ancient state and make it seem a recent relation when it's not.

I can answer that ... I'm afraid no. What you are asking is whether there is any similarity in their cellular receptor binding or "fusion loops" because the motifs look a bit similar.
These are very different viruses and you need to appreciate just how different they are,
Human herpes is a double stranded DNA virus resulting in a chronic infection. It has a vast genome of >150kb, which includes cytokine mimics and proteins that we have no idea what they do. The protein you are referring to is likely a surface antigen in a virus covered in different surface antigens. However, this protein is involved in blocking the beta inferon response is considered a protease - not fusion loop and appears to have some sort of complex replication function involving the capsid. Finally, there is also the question of which human herpes 1 to 8. Here you are referring to HHV1 but for example HHV3 is chickenpox (okay its an acute not chronic infection) and HHV5 is cytomegliovirus (congenitcal infection).
COVID-19 is a simple virus of a 30Kb genome of a single stranded positive sense RNA, causing an acute infection. I agree the Spike protein is much bigger than most RNA viruses, but the receptor-binding is tightly defined with ACE-2. It forms a classic RNA icosahedral virus, singly covered in Spike protein. I have no doubt it will go through conformational change like other RNA viruses in receptor binding but the protein involved is singularly involved in cell entry - thats its job.
Summary
There are large differences between HHV1 and COVID-19 for,

Viruses; DNA double stranded vs RNA single stranded
Infections; chronic vs. acute
Transmission; vertical/integration? direct contact vs. respiratory
Ultimately very different proteins
Finally very different functions

At best this would be heavily  scrutinised as a virological assumption.
Bioformatics If you really wanted to follow this up you could compare the structures around these motifs, because COVID-19 is close to SARS and assess this possibility ... but you would need to be careful about the conformational changes Spike goes through in receptor-binding.
Advice before investing large amounts of time in a bioinformatics question, it is well-worth understanding the biological rationale for the question. In this instance that rationale appears to be lacking. If you found a comparison for example with gamma-coronaviruses, which was absent from many beta-coronaviruses that would be interesting. Here there is simply too much possibility of a chance similarity and given the size HHV1 that cannot be easily overlooked. However it is good that you are thinking and questioning, if you apply the same techniques, e.g. psi-blast, protein-blast into different questions ... that can work.


text:
Identification of starting point for the synthesis of a protein
identify where is the starting point for the synthesis of a protein in the DNA
For indentify a starting point in RNA or DNA my function search all ATG or M in the sequece and report the linear position of start, than the function search TAA and TAG for stop codon and report the end position.
For the end codon what I have to search? only TAA and TAG or TGA and TGG too?
How is possible some protein begin with "F" as in "FFRE..." and not "M..." ?
Than it result more starting point solutions respect to protein found in GeoBank, what are the add restriction for limit that number?
Are there easy to search sequences that make sure that is a start of the protein code? as seems to me "AGG AGG" but is not clear where they should be in the array and not seems general seen proteins.
Can the synthesis of protein be affected by temperature? What other factors come into play?
Passing my function to DNA of one virus, the first test, the function identifies
in the end of virus many little proteins of 2 3 4 5 10 protein letters. These little protein pieces are really sintetized? or there is a min lenght for sintetized proteins?

This seems like a homework question. If you want a software to identify ORFs for another project, @PPK's answer is a good way to start. If you need to make a function to do so as part of a project for a subject, you should look for all stop codons as possibilities.

You seem to be interested in identifying potenital coding regions in a DNA sequence. These regions are called open reading frames (ORF) because they are open to translation.
As detailed in the linked Wikipedia article, an ORF usually starts with an ATG and ends with one of the stop codons (TAA, TGA or TAG in DNA space). However, there are many ORFs that are not translated at all, particularly short ones.
Since finding ORFs is an old and established task in bioinformatics there are many algorithms availabe to do this. I would suggest using such an existing tool because they often have additional features and bugfixes. There is a list in the link above.
Estimating the stability and likelihood of a particular amino acid sequence is much more complicated. There are many tools for protein structure prediction that may be worth having a look at.


text:
How to fit all genes (labels) in chromosome ideogram plot made by RCircos package?
I am using Rcircos to make a chromosome ideogram plot for my gene list (n = 45). However, I am getting this error:
 Not all labels will be plotted.
 Type RCircos.Get.Gene.Name.Plot.Parameters()
 to see the number of labels for each chromosome.

I am using the below-given script and what to know is there any way to fit all my genes in the chromosome ideogram plot
 library(RCircos)
 data(UCSC.HG19.Human.CytoBandIdeogram)
 RCircos.Set.Core.Components(
 cyto.info=UCSC.HG19.Human.CytoBandIdeogram,
 chr.exclude = NULL,
 tracks.inside = 10, tracks.outside = 0)

 #Making a plot with RCircos
 out.file <- "RCircos_IF.pdf";
 pdf(file=out.file, height=8, width=8, compress=TRUE);
 RCircos.Set.Plot.Area();
 par(mai=c(0.25, 0.25, 0.25, 0.25));
 plot.new();
 plot.window(c(-2.5,2.5), c(-2.5, 2.5));

 # plot chromosome ideogram
 RCircos.Set.Plot.Area();
 RCircos.Chromosome.Ideogram.Plot()
 gene.lable.data = read.table ("IF.txt", sep="\t", header=T, row.names = 1) # IF.txt is my input file
 RCircos.Gene.Connector.Plot(gene.lable.data, track.num = 1, side = "in");
 track.num <- 2;
 RCircos.Gene.Name.Plot(gene.lable.data,name.col = 4, track.num)

 dev.off()

My input file looks like
    Chromosome  chromStart  chromEnd    Gene.name
    1   chr1    150363091   150476566   RPRD2
    2   chr1    150549369   150560937   ADAMTSL4
    3   chr1    91949371    92014426    BRDT
    4   chr1    31365625    31376850    FABP3
    5   chr1    150960583   150975004   CERS2

To change the maximum number of genes in RCircos, we could modified the char.width according to this link.
The actual value in a RCircos session could be modified with get and reset methods for plot parameters.
params <- RCircos.Get.Plot.Parameters() #$char.width
params$char.width <- 100 #default 500
RCircos.Reset.Plot.Parameters(params)

#the maxLabels are updated accordingly
RCircos.Get.Gene.Name.Plot.Parameters()

Just want to add that apart from params$char.width, we can also use
params$text.size
params$text.size <- 0.2 #default 0.4
params$base.per.unit < 50000 #imp #default value = 30000


text:
How to add cluster name to metadata in Seurat?
I'm working on a Seurat object and want to name the clusters according to 2 values alone (yes/no). So I want to add a new column to metadata and annotate the clusters (UMAP) with it.
head(B@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito RNA_snn_res.0.6 seurat_clusters
BC01_02       BC01   999789.2         9225 0.2385090196               6               6
BC01_03       BC01   999776.9         9568 0.2925109851               6               6
BC01_04       BC01   999595.1         7482 0.0009888404               6               6
BC01_05       BC01   999992.0         8301 0.3801925053               6               6

I want to add a new column imported from a csv file that contains the row.names same as the meta.data with the status column as yes and no.
head(status)
               status
BC01_02           yes
BC01_03           yes
BC01_04           yes
BC01_05           yes

After adding the column, I want to name the clusters with yes or no, how is it possible?

To add the column status can try the following:
B@meta.data$status <- status[row.names(B@meta.data),]$status

To rename the clusters you can try the following:
B@meta.data$seurat_clusters <- B@meta.data$status

However, I think the plotting functions in Seurat do not use the cluster information from B@meta.data, rather they use information from B@active.idents.
So, I guess you have to do the following instead:
B@active.idents <-  B@meta.data$status

This should work, but, I will verify and update this answer later.

Thank you! It worked!  B@active.idents worked for renaming the clusters!


text:
How to calculate the number of peaks that are upstream/downstream of some other peaks
I have 3 histone marks,I have used Macs2 for peak calling and diffBind to analyze the peaks. I was wondering if you know any way to calculate the peak numbers of one specific histone mark that occur for example within 2kb of peaks belong to another histone mark.
Much appreciated

You can use windowBed from bedtools with -c if you just want counts. Or use -wa -wb and then use groupBy for any other operations.


text:
How to subset genes and its nested features from a GFF file using a gene list
I would like to subset a GFF file (gene and nested features) from a gene list.
The GFF file looks like this
##gff-version 3
Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903

And a have the target genes in a list such as
Genecv11033552m
Genecv11003131m
Genecv11036683m
Genecv11012576m
Genecv11003654m
Genecv11012587m

I know that is possible to subset the gff using grep -f gene_list.txt <gff_file>. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five_prime_UTR, CDS, three_prime_UTR).
This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.
Any ideas?

I think this simple script that uses pyranges solves it. Half of the below is just setup to make a minimal reproducible example.
# pip install pyranges
# or
# conda install -c bioconda pyranges

import pandas as pd
from io import StringIO
import pyranges as pr

#########
# setup #
#########

gff_name = "gtfo.gtf"
contents = StringIO("""Scaffold_1   JGI gene    22901   45904   .   +   .   ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 24031   24062   .   +   2   ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1  JGI gene    22901   45904   .   +   .   ID=Gene2.g;Name=Gene2.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ""")

pd.read_table(contents, sep="\s+", header=None).to_csv(gff_name, sep="\t", index=False, header=False)

genes_to_keep = set(["Genecv11000001m.g"])

############
# solution #
############

df = pr.read_gff3(gff_name, as_df=True)

def all_descendants(df, genes_to_keep):
    old_len, new_len = -1, len(genes_to_keep)
    while old_len != new_len:
        genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
        old_len, new_len = new_len, len(genes_to_keep)
    return df[df.ID.isin(genes_to_keep)]

desired_gff = all_descendants(df, genes_to_keep)
print(desired_gff)
#    Chromosome Source         Feature  Start    End Score Strand Frame                                ID               Name longest             Parent
# 0  Scaffold_1    JGI            gene  22900  45904     .      +     .                 Genecv11000001m.g  Genecv11000001m.g     NaN                NaN
# 1  Scaffold_1    JGI            mRNA  22900  45904     .      +     .                   PAC4GC:50510902    Genecv11000001m       1  Genecv11000001m.g
# 2  Scaffold_1    JGI  five_prime_UTR  22900  23284     .      +     .  PAC4GC:50510902.five_prime_UTR.1                NaN     NaN    PAC4GC:50510902
# 3  Scaffold_1    JGI             CDS  23284  23423     .      +     0             PAC4GC:50510902.CDS.1                NaN     NaN    PAC4GC:50510902
# 4  Scaffold_1    JGI             CDS  24030  24062     .      +     2             PAC4GC:50510902.CDS.2                NaN     NaN    PAC4GC:50510902

This is super nice, and I believe way simpler than the script I had come up with! Many thanks!

My scipt outsources the parsing of the gtf to a library, which is probably why it looks nicer. Much complexity is hidden behind that read_gtf call XD

Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython.
Dealing with nested data in bash or with some fancy awk script is probably going to trickier to write than it's worth.
You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with mRNA in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.
If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.
Repeat until you reach the end of the GFF file or your list of target genes is empty.

I have actually done it in python using the "flow" you suggested. Thank you!

Answer as a python script for this task specifically

That script looks good! My only suggestion is that, depending on the GFF file you're loading, you may get slowed down when trying to load the entire file as a pandas DataFrame. If you are trying to parse a multi-gigabyte GFF file, going line-by-line keeps the memory usage low. Otherwise, what you've done looks good

Precisely, I have already felt that pain :). I am actually working on a better solution for this. Thanks for the help!


text:
Standard for when to collapse phylogenetically uninformative nodes
I have been creating phylogenetic trees using RAxML (using Boostrap + ML) and I end up with a bestTree file, outputted in newick format. When I open up the tree in FigTree, I notice that it collapses many of the internal nodes with small values (1 x 10-6). I am working on visualizing the trees using ETE, which when I render the tree includes these small breakpoints. I was wondering if there is a standard for when to collapse nodes that are considered phylogenetically uninformative that I can apply. Thank you for any help.

Hi @Ace8910 what do you mean by E-6 .. are you referring to genetic distance?

If 1E6 is genetic distance, this means one nucleotide difference per million base pairs, which is of course incredibly small and bootstrapping at this level doesn't really work. So what I suspect is happening is that the bootstrap is so low it is not present on the tree. You can assign a cut-off bootstrap but anything below 75% is not interpretable in any case. What might be happening is you are producing the bootstrap "cladogram" which is showing a default branch length, but when <1E6 genetic distance it is showing a trifurcation. This simply means the bootstrap of the consensus tree is so low it is collapsing the branch and this is a function of the very small genetic distance at this node.

Looking at your treefile I don't see any collapsed clades.
(Group_9:0.00000100000050002909,(Group_50:0.00507659947127924935, (.

The clade above has not been colapsed it is just a very small distance in a paraphyletic group. It looks fine to me, i.e. Group_9 is 1E6. It could look 'collapsed' simply because it is so small it is difficult to see. I hope thats ok?

Could you post post of the raw tree file to examine this in more detail?

My apologies, it is collapsing nodes where the genetic distance to that node is 1e-6. Here is a snippet: (((Group_40:0.00000100000050002909,(((Group_14:0.00506368298581611701,Group_5:0.00504476117481173891):0.00000100000050002909,(((((Group_6:0.00506368450410310620,Group_18:0.00503828136153215873):0.00000100000050002909,(Group_1:0.00506367795420949521,(Group_7:0.00000100000050002909,(Group_9:0.00000100000050002909,(Group_50:0.00507659947127924935, (..........       I can also provide the tree file by email if that would be helpful.

It looks fine from the treefile you posted (above)

Thanks Michael, yeah the sequences were only 500 nt long, so I believe those breakpoints aren't informative in terms of mutations occurring. We are trying to map instances of mutations to a reference genome using the trees created.

Thanks, please tick my answer as "accepted"  if it was useful. We can deal with other other questions separately on another day. You may need ETE3 when you scale up to automate the your search.


text:
How identifiable are human omics data and how to mitigate their identifying features?
Say a database were to store human omics datasets. The human subjects are known and the sample size is rather small in size initially (n=500). The database contains genomics, transcriptomics, proteomics, gut metabolomics, and epigenomics. Since the sample size is rather small initially, it would be important that certain identifying features are mitigated so that individual subjects cannot be identified. This leads me to two questions:

What identifiable features could be present in these different types of raw data? (Such as race, sex, age, hair color, eye color, height, where someone has lived, and anything else I may not even be considering)? Which of these omics types would be the most and least dangerous for identifying individuals? How can the identifiable nature of these data be mitigated?

If environmental metagenomics is also collected by these human subjects, is it possible to identify the human subjects by contamination? (i.e. some of the reads from the metagenomics data inadvertently contain human reads?) How can the identifiable nature of these data be mitigated?

I think this subject may be a bit futuristic, but I am very uninformed. If there are any references that provide additional thinking about these topics, please kindly share. Thank you for sharing your thoughts.

Its not futuristic, it has almost certainly been done both the db and the metagenomics. However @clove444 the site works best by asking a specific technical question.

This is one of the major problems with genomic information in todays research. This was highlighted some years ago with the police using publicly available genomic data bases to idenitfy unkown murder suspects.
The full extent of this issue was exemplified in a
paper by Yaniv Erlich from 2018 that should be a good starting point for you.
They claim that giving the current (2018) amount genetic information more than half the searches for a person of european descent will result in a thrid cousin or better match. And these odds will only inprove with time.
Therefore, the genetic information not only can give you a hint at inherited variables like the approximate height of a person but can even yield their name and address.
About your second question, I would think that the amount of host (i.e human) material collected alongside the microbiome sample is rather small. This probably is not sufficient to obtain enough coverage for a profile. But I could be mistaken.

PPK provides a great answer, but for question 2 I can provide a different perspective. For shotgun metagenome sequencing (without any enrichment/depletion protocols) it is common for >90% of reads to map to the human genome. There is variation across body sites, for example the gut microbiome is rich and so the % of human reads from a stool sample will be less compared to say a blood sample.
Checkout this study, where they mined human microbiome datasets for host (human) reads and were able to re-construct draft full human genomes with sometimes 20X coverage. See also part of their discussion:

We show here that it is possible to reconstruct complete host genomes
using metagenomic sequence data, which is potentially identifiable.
However, this was possible due to the unique study design of the HMP,
whereby multiple body sites from each individual were sequenced at a
high depth, allowing us to pool data across body sites and reach a 10x
mean coverage per host genome. Common metagenomic shotgun sequencing
studies, which usually include an order of magnitude less sequence
data, are unlikely to enable such an analysis. Moreover, the majority
of studies sequence stool samples, which include many fewer
host-derived reads. Nevertheless, we anticipate that future shotgun
metagenomics sequencing studies would consider these potential privacy
concerns.