text:
Tool for rna/lna melting temperature prediction
Is there any tool available to predict the melting temperature for rna/lna oligos with its perfect complement ?
I know this amazing website from qiagen/geneglobe that estimates the tm value.
But is there any command line tool that does this?
I can estimate the normal tm value with RNAcofold and the -T option but there is no lna option.
I searched for a while but i can only find online services.

There are several programs for simulating/predicting DNA melting but you could try them for RNA.
You can look here http://bioinformatics.org/meltsim/wiki/Main/HomePage.

Both Bio.SeqUtils.MeltingTemp of Biopython and Primer3 offers different methods/algorithms to calculate the Tm value.


text:
choosing a good representative genome subset
I'm trying to build a genomic database for DNA alignments.
I started with NCBI accessions, but the data is very multiplicative, so I want to use subset of [max] N different strains for each specie.
my question is: I can i obtain the minimal N that the samples will be diverse enough and will cover the specie good?

Depends on what you are looking at (species) and secondly what post-alignment analysis you are planning.

Check out the progenomes database: http://progenomes.embl.de/

source: Stack Exchange

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

sort: votes


text:
Wrong blast results
I'm sorry if this is a naive question, but I get an odd thing going on.
I'm performing blast remotely in this way:
'blastn -task megablast -outfmt "6 std ssciname" -reward 1 -penalty -3  -db nt_v5 -max_target_seqs 100 -word_size 22 -query ' + final_file + " -taxidlist " + Plants_DB + " -out " + blastfile

However, when I see my remote results, I get a species whose seqid gives me Vaccinium oxycoccos. And when I put this seqid in blast online, the species is totally different.
(With this blast command, I perform a blast search against a specific list of IDS.)
How can this happen?
Thanks

Are you sure the parameters including word size and database, as well as all default settings are identical?

What exactly do you mean by seqid? Species have taxid, e.g. humans are 9606. Whereas genbank files have accession numbers (and formerly gi numbers).

@MatteoFerla I have this result, 4F703:00578:02013 KX677575.1 100.000 320 0 0 1 320 700 381 3.25e-165 592 Vaccinium oxycoccos and when I put the second result in blast I get a different specie...

@MichaeG. I tried to specify the most part with this commands: blastn -task megablast -word_size 28 -reward 1 -penalty -2 -outfmt 6 std ssciname -db nt_v5 -max_target_seqs 100 -query input.fasta -taxidlist test_ID_plants.ids -out blastfile.txt... do you think there is a problem with my nt database?

what do you mean with remotely? because in your command you are not using the --remote option, and there is an extra single quotation mark as well.

@zorbax, you don't need the option remote necessarily. also, what do you mean by "extra single quotation mark"

I assume that you call "remotely" because you are running your process in a remote computer. But the -remote option in blast means remote service at the NCBI servers, otherwise you are using a local version of nt database called nt_v5.

If I add the remote option to my inital command, the result continues to be wrong...

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


text:
Given a multiple sequence alignment, how do I output all of the non-consensus characters by location?
I have an MSA from MAFFT (as a FASTA file; there are several sequences in the MSA), and I would like to quantify the "variants" between these sequences.
What is the standard tool to do something like this?
I would prefer to know the non-consensus variants in the the entire MSA, but that may be too tricky, so I would like to figure out the non-consensus variants between those select sequences.
I'm familiar with calculating an MSA, and then visualizing it---but I would like an output file with all of the different bases.

This is surprisingly not done or requested in the community, I'm also struggling.

A very basic solution in python:
from Bio import SeqIO
seq_dict = {i : rec.seq for i,rec in enumerate(SeqIO.parse("your_msa.aln", "fasta"))}
for pos in range(len(seq_dict[0])):
    pos_bases = ([seqio[i] for seqio in [*seq_dict.values()]])
    # If position contains more than one unique base
    if len([*set(pos_bases)])>1:
        print(pos, pos_bases)        # prints the position and alternative bases there

If you need something more useful, just leave a comment here and I'll develop this further.


text:
Why index a sequence file with minimap2?
I have been using minimap2 for a few month now and I just found a script that follows these steps:
minimap2 [-x preset] -d target.mmi target.fa
minimap2 -a [-x preset] target.mmi query.fa > output.sam

I have never seen this .mmi index before. I always used minimap2 directly like this:
minimap2 -a [-x preset] target.fa query.fastq > output.sam

So my question is why should I index my reference genome before mapping? Is it only for speed up? Will it up speed my code if I do it? Or it is already done without specifying it?

If you don't provide an index, minimap2 will just generate the index on the fly (so, yes, it's done without specifying it). As such the mapping step will take longer. However, index generation is supposed to be very fast in minimap2 (a few minutes for the human genome), so it probably won't matter much.
Doing it separately will therefore also not speed up your code unless you end up using the pre-generated index multiple times.


text:
Optimum parameters for genome assembly via canu
I have been trying to assemble my genome via canu. Canu gave me about 200 contigs using default parameters. My genomesize is ~9Mbp. Can you recommend some parameters I can tweak to decrease the contigs number?

You can also try other de novo assembly tools like SMARTdenovo to see if they output results more in line with what you expect

The canu parameter reference guide is always a good place to start.
It would help to know a little bit more about what your assembly looks like. E.g.

what is your n50, n90...?
What is your dataset coverage?
Have you run any tools like quast to evaluate quality?
What is your input data type (probably PacBio)? Read length distribution?
What is your expected chromosome/replicon number?
Have you compared the assembly to any other data, e.g. genomes of close relatives? (quast can do this)

Basically, why do you think that 200 contigs is too many contigs? Sometimes you just don't have enough coverage to get more contiguous/"nice" assemblies. That's not to say it's not possible, of course, just that we don't have enough information to evaluate.
It's possible that some of those contigs are small and artefactual, so you might get rid of them with just length/coverage filtering. You can do some of this in canu itself (see "Output Filtering" section).

It is pacbio reads. It results in a contig which is circular and ~9.5 Mbp. Other contigs are either 2Mbp or lesser. They do not seem to be plasmids as well from annotations, and their coverage is also good, ~20X (for the 9.5Mbp contig, its ~30X). Also, other contigs are not repeats of the 9.5 Mbp contig. The longest contig is similar to the reference genome from mauve alignment and dot plot.
I am confused if even the 2Mbp contigs could be artefactual?

@menten well, to me it sounds like your circular contig is your reference genome, congratulations. the fact that you get other contigs in addition to it suggests to me 1) you have contamination. 2) the other contigs are haplotigs (contigs from alternate haplotypes or clones in your culture) of your genome. I would suggest doing an all-by-all alignment of your contigs via minimap2, which should identify any haplotigs. do you know the ploidy? if it's diploid and not homozygous, well, then you likely have haplotigs with high-coverage pacbio.

(i am guessing that you are working on bacteria in which case I know high ploidy may be an odd question, but you never know)


text:
a package with randomization/permutation tests for genomic analyses
please would you let me know, is there a package in R/python/C that has implemented randomization/permutation tests for joint genomic analyses
(an example of joint genomic analyses -- when jointly considering both GENE EXPRESSION and PROTEIN BINDING along the DNA).
The context of my question is the following :
let's consider 1000 UP_REGULATED genes with increased PROTEIN X and with a HISTONE MARK Y :
in order to show that PROTEIN X is related to HISTONE MARK Y for 1000 UP-REGULATED genes, what "controls" would you use for the comparison :
-- 1000 RANDOM GENES (and multiple randomization tests)
-- 1000 UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y
-- 1000 UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y
-- 1000 UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y
-- anything else ?
any suggestions are very welcome. thanks a lot,
bogdan
ps : 've posted a similar question on BioC website : https://support.bioconductor.org/p/128828/

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

https://support.bioconductor.org/p/128828/

I would simply look at co-variance across the entire dataset and examine values >0.8 in the first instance rather than select 1000 genes.
With your approach sampling 1000 genes can be criticised for selective sampling and the way to get around this is to perform replicates, with (bootstraps) or without replacement. You could produce e.g. 1000 replicates of 1000 genes and perform co-variance on each replicate to examine Y = X.
Python has a specific permutation/combinations module call itertools, which is different to random resampling and doesn't need replicates
from itertools import combinations
X = [1,2,3]
Y = [1,2,3]
print (list(combinations(iterable, 1000) + X + Y))

The standard way is to sample is,
from random import sample
mygenes = [1, 2, 3, 4, 5]

print(sample(mygenes,3))

Output
[2, 4, 5]

There is nuance and further funcationality ... I cannot remember is "sample" is with or without replacement (its quite a big deal)

thanks a lot, Mike ! very helpful ! wondering if possible to recommend any package in R has the functions for bootstrapping, or similar functions as in itertools. thanks a lot !

Hi @bogdan, I don't use R so can't recommend anything.


text:
Running MaxQuant on Linux
I'm trying to run MaxQuant on a linux laptop, hower I'm constantly running into problems with MaxQuant crashing at different stages
I tried to use

mono 6.8
mono 5.4.1

on

centos
ubuntu

each with MaxQuant

1.6.3
1.6.7
1.6.11

(everything using docker), but I had no success so far.
I also tried to use the ahead of time compilation of mono, also without sucess.
MQ mostly stalls at the "preparing searches" stage.
Has anyone managed to get MQ running properly on linux and can tell me, which linux distro/version, mono version and MQ version you are using?
Any other help is of course also welcome

There is a conda package for maxquant that should work on most platforms (I've used it on CentOS and I think it's been used on Ubuntu as well). There's also a docker container version of that that we (Bioconda) have posted to the biocontainers project that should also work.
Please note that large parts of a maxquant analysis are single threaded, regardless of what you specify. We've also had single analysis runs lasting over a month (even when using 90 threads), so the program can be very slow if your search space is large.

I have been running MaxQuant quite frequently on Linux. I would use conda to create a dedicated environment and install mono.
# create the environment
conda create -n maxquant -c conda-forge mono
# activate the environment
conda activate maxquant
# run any maxquant version
mono /path/to/maxuant/MaxQuantCmd.exe mqpar.xml

This installs a conda environment with mono which is a .NET framework for Linux. I had problems in the past running certain MaxQuant versions that way.
The mqpar.xml that is required can be created with the MaxQuant GUI. And then the path for input files needs to be contained correctly. The path to the fasta file and the raw file. I would create a template file and then parse these filenames into the specific locations.
# mqpar.xml
...
<fastaFiles>
  <FastaFileInfo>
     <fastaFilePath>__FASTA__</fastaFilePath>
     <identifierParseRule>>([^\s]*)</identifierParseRule>
     <descriptionParseRule>>(.*)</descriptionParseRule>
     <taxonomyParseRule></taxonomyParseRule>
     <variationParseRule></variationParseRule>
     <modificationParseRule></modificationParseRule>
     <taxonomyId></taxonomyId>
  </FastaFileInfo>
</fastaFiles>
...
<filePaths>
  <string>__RAW__</string>
</filePaths>
...

The combination of maxquant=1.6.10.43 and mono=5.14 worked and a few other versions. Some versions between 1.6 and 2.x did not work at all.
Specific versions of MaxQuant are also available directly via conda and there is no need to download the zip file. Instead, you can simply install it with:
conda install -c conda-forge -c bioconda maxquant mono=5.14

Recently MaxQuant 2 has been added, and not all previous versions are available and run maxquant with:
maxquant mqpar.xml


text:
Pol_II genome wide occupancy analysis
I would like to know how to perform this global Pol II occupancy which is shown here:
 Paper
Figure 1G
So in this figure they are analysing the Pol II occupancy in Control,PAF1 KD and LEO1.
I have similar data.

In order to determine whether PAF1 has a direct effect on the regulation of Pol II at the promoter-proximal region, ChIP-seq was performed. Heatmaps ranked by decreasing Pol II occupancy show a positive correlation with PAF1 and Pol II (Figure 1G). To corroborate these findings, we also performed ChIP-seq with LEO1, as this was the only other antibody that we tested that gave any significant ChIP-seq signal. Like PAF1, LEO1 is found throughout gene bodies with particular enrichment around the TSS (Figure 1G). Interestingly, LEO1 was not broadly lost from chromatin upon PAF1 knockdown, but the distribution was altered, including reduced occupancy in the promoter-proximal region

So what I have done so far in terms of analysis is I did peak-calling using homer for my control and KD condition along with their Pol II chip.I was going for this differential binding analysis but this is now how i Pol II data is done after reading the paper I got. In homer there is this function analyzeRepeats.pl which is provided in this section Measuring Gene Expression in Exons vs. Gene Bodies. telling genes (default) - Counts tags on the full gene body (TSS to TTS).This is useful for GRO-Seq where we expect coverage across the entire transcript.Can also be used to quantify H3K36me3 or PolII ChIP-Seq.
How do I proceed and what tool should I use to show the Pol II occupancy at different condition and as of now I don't have respective rna seq data, do i need expression data of genes to perform occupancy analysis.

comment discussion between @krushnach-chandra and @stupidwolf (in chronical order):
SW: No you don't need gene expression data.. You can just count the reads from polII chip (normalized by some method of your choice) and perform the analysis
KC: are you suggesting I can compare the Pol II chip between control and Knock out with their respect IP data?
SW: more or less yes...if you are interested in genomewide effects.. For example, if the treatment causes the polymerase to stall, you will see more reads in the promoter compared to the gene body
KC: so naive question do i need those Input and IgG at all for the comparison?
SW: good question.. depends on how you analyze it. If you always normalize within the gene, for example, you calculate a coverage normalized to that in the gene, or some meta gene thing, you don't need.
SW: If you do a differential gene expression based analysis, in theory you don't really need that, because you will always compare between your two conditions right. The IgG is good for checking you don't have any odd stuff... You can see whether there are genes where the control is enriched over your samples
KC: "calculate a coverage normalized to that in the gene, or some meta gene thing" any example can you cite me..so that i can follow and proceed ,as it looks like binding data are really complicated compared to rna seq
KC: "If you do a differential gene expression based analysis," you saying i can take out counts and do diff analysis for chip seq data using deseq2? but how do i get the information regarding if Pol II activity is going or stalled ? and how do i find the occupancy from the differential result
SW: welcome to the world of epigenomics!!! maybe something like this
SW: As for differential gene expression, yes. If you do this separately for promoters and gene body, you can see that the fold change will be different for a gene's promoter vs gene body, in situations where there's stalling, or reduced elongation etc.. I think it's a good way to account for the variation
KC: your insights are good start for me to run the analysis
KC: mundane question how do i get promoter and gene body information is it after annotation which should be done after running analysis ?
SW: Meaning? I normally align the reads, then for TSS and gene body you have to define it yourself using the transcript or gene annotation. Normally like +/- 200 bp around TSS for promoter and gene body is rest of it?
KC: "Normally like +/- 200 bp around TSS for promoter and gene body is rest of it" so i can do that in both chipseeker as well as homer?


text:
Parsing a GenBank file with multiple gene entries
Basically a GenBank file consists of gene entries (announced by 'gene') followed by its corresponding 'CDS' entry (only one per gene) like the two shown here below. I would like to extract part of the data from the input file shown below according to the following rules and print it in the terminal.

There are two blocks of gene data shown below. I want to extract part of both blocks.
The extracted text for each block starts with a line that contains spaces at the beginning of the line followed by gene
The extracted text for each block ends with a line that contains /db_xref="GeneID
The perl and awk tags are just suggestions. An answer can use a different program(s).

Input file (source: biotech):
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
                     /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"
                     /translation="MQTVNVDVAIVGAGGGGLRAAIAAAEANPNLKIALISKVYPMRS
                     HTVAAEGGAAAVAKEEDSYDKHFHDTVAGGDWLCEQDVVEYFVEHSPVEMTQLERWGC
                     PWSRKADGDVNVRRFGGMKIERTWFAADKTGFHLLHTLFQTSIKYPQIIRFDEHFVVD
                     ILVDDGQVRGCVAMNMMEGTFVQINANAVVIATGGGCRAYRFNTNGGIVTGDGLSMAY
                     RHGVPLRDMEFVQYHPTGLPNTGILMTEGCRGEGGILVNKDGYRYLQDYGLGPETPVG
                     KPENKYMELGPRDKVSQAFWQEWRKGNTLKTAKGVDVVHLDLRHLGEKYLHERLPFIC
                     ELAQAYEGVDPAKAPIPVRPVVHYTMGGIEVDQHAETCIKGLFAVGECASSGLHGANR
                     LGSNSLAELVVFGKVAGEMAAKRAVEATARNQAVIDAQAKDVLERVYALARQEGEESW
                     SQIRNEMGDSMEEGCGIYRTQESMEKTVAKIAELKERYKRIKVKDSSSVFNTDLLYKI
                     ELGYILDVAQSISSSAVERKESRGAHQRLDYVERDDVNYLKHTLAFYNADGTPTIKYS
                     DVKITKSQPAKRVYGAEAEAQEAAAKKE"

Desired output:
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

I tried using pcregrep --multiline .*'START-SEARCH-TERM.*(\n|.)*END-SEARCH-TERM' path/to/SOURCE-FILE. I installed pcregrep (grep utility that uses Perl-style regexps) in Ubuntu with sudo apt install pcregrep. I also installed Biopython with sudo apt install python3-biopython and ran the Simple GenBank parsing example from Biopython Tutorial and Cookbook.

I would strongly suggest simply using biopython, bioruby or biojulia etc. These libraries are really good for extracting data from genbank files. You tagged perl, bioperl can definitely do it, but Perl is a dying language and bioperl is not straightforward to install IIRC from 10 years ago —but worth a try.

@MatteoFerla take that back! Just because young whippersnappers today don't appreciate the power and beauty of Perl does not make it a dying language! Them's fighting words! :P

Yeah agreed, code is code. Two things will continue Perl in any age, regex and Perl one liners (definitely stylish). Have you ever heard of a Python one-lliner?

You can simply use grep for this purpose as shown below,
grep -v /translation bio.txt | grep -B100000000 /db_xref= > output_file.txt

Just make sure that you keep the number with B bigger than the number of lines of your file.
If you print the contents of the above file you get your desired output as given below,
:~$ cat output_file.txt
gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

You can do this very easily with awk:
$ awk '/^ {5}\w/{a=1} /\/translation/{a=0}a' file.gb
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

The idea here is to set a to 1 if this line starts with 5 spaces followed by a word character. Then, we set a back to 0 if this line matches /translation. The default action for awk when an expression evaluates to true (not 0) is to print, therefore the final a will cause all lines read while a is not 0 to be printed, effectively removing everything after each /translation line.
A more easily understandable version of the same code would be:
awk '{
        if(/^ {5}\w/){
            a=1
        }
        if(/\/translation/){
            a=0
        }
        if(a==1){
            print
        }
    }' file.gb

@karel it prints /db_xref="GI:219870283" but not /db_xref="GeneID:7278620".

@karel see update. Awk to the rescue!


text:
Extracting part of a csv by a list
I have a .csv file like this
pair,S17A,S2..APOBEC.,S1..age.,S18.like,S3..BRCA.,S17B
LP2000104-DNA_A01_vs_LP2000101-DNA_A01,0.448027080410302,0,0,0,0,0.410602648202438
LP2000325-DNA_A01_vs_LP2000317-DNA_A01,0.343096214261112,0,0.111190317814775,0.171937435363758,0.187395051556141,0.156209735124657
LP2000327-DNA_A01_vs_LP2000319-DNA_A01,0.436249785689912,0,0.066118729867486,0.077392575088635,0.116384601626772,0.263251836199875

I have a list of samples which are a part of first column for example
LP6005334-DNA_H01
LP6008334-DNA_A03
LP6008334-DNA_A04
LP6008334-DNA_B02
LP6008334-DNA_C02
LP6008334-DNA_D02
LP6008336-DNA_H01

I want to extract .csv for this list of samples
I can not think of intersect function as these samples are a part of sample name not full name
Can you help?

As usual you do not provide example output despite in like every of your question people ask for it. No surprise you constantly get downvoted.

You can use strsplit.
ciao <- "LP2000104-DNA_A01_vs_LP2000101-DNA_A01"
sampleIneed <- "LP6008334-DNA_A03"
sampleIneed %in% strsplit(ciao,"_vs_")[[1]]
[1] FALSE

"LP2000104-DNA_A01" %in% strsplit(ciao,"_vs_")[[1]]
[1] TRUE

You can loop (e.g. sapplyor lapply) through the element of the first column to extract only the rows you need.
You can also simply use grep on the first column
ciao <_ your_csv
grep("LP2000104-DNA_A01", ciao$pair)


text:
Download sequences of isoforms
I want to collect all isoforms of all genes (Fasta/Fastq file of nucleotide and protein sequences) Arabidopsis Col-0. I am wondering if there is a straightforward way to download the file from any database.
I am new to sequence analysis, so I apologize if asking a very basic question.

Do you want a single gene? Look at Uniprot then

IMO you are into alternate splicing sites, more proteomics and RNAseq than genomics.

Col-0 is, to my knowledge, the standard reference Arabidopsis thaliana ecotype, so the easiest way to get the sequences for all genes is to just download them from Ensembl. The nucleotide sequences are available here and the protein sequences here.

For RefSeq data, go to NCBI homepage and search for Arabidopsis Col-0 genome. Click on the 'Get Data' button in the result box and you will be directed to the NCBI Datasets portal with the latest RefSeq assembly selected. Click on the blue 'Download' button, uncheck the genome and GFF3 files from list, check transcript and protein FASTA files from the list and download data.


text:
Merge multiple PDB files into a single file
I have several pdb files and each of them corresponds to a domain on a protein. I want to merge these files into one single pdb for the whole protein, and the approach I'm looking for should allow me to order these domains according to the order by which they're located on the protein (for example, merging domain A, B, and C by the order CBA). It would be preferable if I can do the merging with VMD but any suggestion would be appreciated.

PyMOL is a rather flexible tool for doing this.
The command create combo, domA or domB or domC in the PyMOL app will create a single object. However, they will need to be the same chain and segment (segi). Also you will need to fix the numbering to match what you expect them to be.
To fix this you need to do:
 alter domA, chain='A'
 alter domA, resv+=0
 alter domB, chain='A'
 alter domA, resv+=100
 alter domC, chain='A'
 alter domA, resv+=200
 alter all, segi=''
 sort
 create combo, domA or domB or domC
 save combo.pdb, combo

PyMOL can be installed as a Python 3 module, which works in the same way but with pymol.cmd and is a lot more replicable. But for starting out the App is fine.
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.alter('domA','chain="A"')
     etc.

Now. This assumed that the coordinates are where we want them to be. Which is highly unlikely and that there are zero overlaps.

If there are overlaps, aligning the overlapping residues (align domB and resi 80-100, domA and resi 80-100) can be done. Just remember that the first element is the mobile object, the second the target, hence why I wrote it backwards.
If do not have overlaps and want the domains docked together you can use the DEMO tool by the I-Tasser people.
If you want your domains as "balls on a string" aligned on the x-axis for a nice picture and you do not care if it may not be biological, you may want to check out a Github repo of mine for doing just so (it offsets, aligns if possible and spaces them out, but does not add linkers).

If you want to add linkers, you can use submit your sequence to Phyre or I-Tasser with your missing density model as a reference, or you can add them with rosetta remodel (tricky) or make them up in Coot or PyMOL (bad idea, but if you are in a hurry... here is a tutorial of mine).


text:
HTSeq-count is returning 0 for every gene, instead of expression value
(Ported from stack overflow)
I'm trying to summarize gene count using htseq-count; and it's returning 0 counts at every gene. I'm not sure what I'm doing wrong; I think it has to do with the gene flag I'm using.
I've tried using the GTF for Arabidopsis TAIR 10: which I got from: ftp://ftp.ensemblgenomes.org/pub/release-46/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gtf.gz as well as the GFF file available on Araport.
I'm using various BAM files I know to have gene expression values.
I've tried both locally and on galaxy
The local commands I've tried so far are:
python3 -m HTSeq.scripts.count --idattr=gene_id #bam #gtf > summary.txt
python3 -m HTSeq.scripts.count -i gene_id #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i locus_tag #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i ID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i SeqID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i Parent #bam #gff > summary.txt

I've done the same on galaxy; inputing the -i value in.
Any help as to where I'm going wrong would be great.
Thanks!

hey @sdshinghal, check a few things. 1. ensure bam file is ok (you can look at it), 2. like the answer below, check your gtf matches the reference genome

Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like 1, some other reference sources use chr1.
Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.


text:
how to generate structure file from blast results of an input sequence
I have tried to incorporate blast in my application using biopython.i could extract the blast results but now i want to extract the structure file corresponding matches. My code is :
from Bio.SeqIO import parse
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

file = open("prats.fasta")

records = parse(file,"fasta")
for record in records:
    print("Id: %s" % record.id)
    print("Name: %s" % record.name)
    print("Description: %s" % record.description)
    print("Annotations: %s" % record.annotations)
    print("Sequence Data: %s" % record.seq)
    print("Sequence Alphabet: %s" % record.seq.alphabet)

from Bio.Blast import NCBIWWW
sequence_data = open("prats.fasta").read()
sequence_data

result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
result_handle

with open('results.xml', 'w') as save_file:
    blast_results = result_handle.read()
    print(blast_results)
    save_file.write(blast_results)

from Bio.Blast import NCBIXML
E_VALUE_THRESH = 1e-20
for record in NCBIXML.parse(open("results.xml")):
    if record.alignments:
        print("\n")
        print("query: %s" % record.query[:100])
        for align in record.alignments:
            for hsp in align.hsps:
                if hsp.expect < E_VALUE_THRESH:
                    print("match: %s " % align.title[:100])

from here i want to generate the xml files and then fetch structure of each nucleotides results

If you the DNA/RNA crystal structure:
Simply change the nt in NCBIWWW.qblast("blastn", "nt", sequence_data) to pdbnt to get the nucleotide sequence in a chain.
The PDB code will be the accession ID with an underscore plus chain letter.
Do note that the PDB structure may be offset.
If you want the crystal structures of the encoded protein:
blastx (translation) on the pdbaa database. As above the PDB code will be as above.
If you want a model of the encoded protein:
Check Swissmodel for a pre-existing close threaded model, else check out I-Tasser (best), Phyre (fast).
If you want a model of the DNA/RNA:
You are running blast without special settings —its not a short motif—, so these things may be long, so do be weary of any result. There are a few RNA folding prediction tools online (e.g. RNAFold).


text:
How to calculate fold change in gene expression from RSEM (RNAseq)
I have gene expression data from RNAseq, specifically: log2(x+1) transformed RSEM normalized counts.
How can I convert these data from multiple samples to determine fold change in gene expression?

You should use a proper statistical framework for RNA-seq dfferential analysis (which includes FC calculation). Standard tools for this are (among others) edgeR or DESeq2. You could use tximport to import RSEM outputs into R and then use its output for e.g. DESeq2. The linked manual provides example code for this. Please read the manuals of edgeR or DESeq2 towards the respective downstream analysis. Both are extensive and compreshensive.

Thanks a lot for this information! So, just to confirm, that means that there is no 'simple formula' that can convert RSEM to fold change? Thanks again.

No. By the way, fold changes without statistics are basically meaningless. Just run it through DESeq2 and then take the results for your analysis.

Thank you very much. I will give this a try.


text:
Practical use of RNA structure
I have spent some time studying RNA structures. While there is a range of interesting phenomena and functions that certainly deserve understanding from the scientific point of view, I have never encountered a practical use of this knowledge (in terms of disease diagnostics, medications, etc.) Have you? If a practical use of RNA structures is still in the future, what could that be?

Hi @Vadim, the Bioinformatics Stack Exchange system doesn't work well for open-ended, or opinion-based questions. If you have a specific problem, please state that instead. Otherwise, I'd recommend asking this question on the Bioinformatics subreddit.

Please look up flavivirus 'double loops' as you described them previously (post for "Coronavirus RNA') and associated RNA secondary structure anomalies for dengue virus and associated vaccine (Butantan) and the yellow fever virus and its vaccine (17D). If you are aware of "double loops", we must be aware of the association of RNA secondardy structure and attenuated strains, i.e. live vaccines. The point is that modification of the RNA secondary structure is a key strategy in flavivirus attenuation.

This is a good point. Is this a widely used technique or is it limited to flaviviruses?

I am not aware of this outside flaviviruses, but my interests are restricted to arboviruses, whilst yours are broader.

I am more on the mathematical/computational side, and I worked mainly on HIV.

When I have a minute I will edit the post to describe the attenuation of flaviviruses via the 3' UTR (untranslated region), this was traditionally about  RNA 2ndary structure modelling of pseudoknots.


text:
What is the aim of insertion codes in the pdb file format?
In PDB file we always see insertion codes (ic) in a few proteins. Could anybody please tell me why it has been added in structure and what is the aim of that ???

Hi @user7355 there isn't enough information to answer this question. An example protein would useful together with information about the protein e.g. resolution (Amstrongs), method of structural determination

@MichaelG. the insertion code has a well defined meaning, it's part of the sequence ID.

When I wondered it myself years ago and I found a good explanation in a mailing list archive. I think it was this one, by Frances C. Bernstein:

The problem with numbering started when people wanted to compare the
  'same' proteins from different species.  They found that there were
  the following possibilities that gave rise to differences:

More or fewer residues at either end.
Extra residues at various places within the chain.
Fewer residues at various places within the chain.
Different amino acids at the same place.

Now imagine that residue PHE 195 is very important for the activity of
  the protein in species A.  But in species B it is residue PHE 197 and
  in species C it is PHE 212, because species B and C are not the same
  length as species A.
Because people felt it was important to preserve the amino acid numbering for 'important' residues and to be able to readily discuss
  and compare the structures from different species, various people
  decided to try to number the proteins from species B, C, etc. to match
  the numbering used for species A.  In doing this, one must have gaps
  (missing numbers) where a sequence is shorter.  But what should one do
  when a sequence is longer?  This is the case where it is necessary to
  insert extra numbers and this is done by using insertion codes.
Thus the insertion code is an integral part of the residue number and it is improper to ignore that field when using a PDB entry.  You
  must also allow for 'missing' numbers when using a PDB entry.

In the PDB format the insertion code is in column 27, following the sequence number (columns 23-26):
                      nnnnI
ATOM  11918  CZ  PHE D 100      -6.852  76.356 -23.289  1.00107.94           C
ATOM  11919  N   ARG D 100A     -9.676  74.726 -19.958  1.00105.71           N
...
ATOM  11970  CE  MET D 100H     -8.264  83.348 -19.494  1.00107.93           C
ATOM  11971  N   ASP D 101     -11.329  81.237 -14.804  1.00107.41           N

Together, columns 23-27 make a sequence ID, which is a mostly-numeric ID of the
residue position.
Note: It's not a full ID of the residue. When you have a point
mutation (a.k.a. microheterogeneity) you have 2+ residues with partial occupancy at the same position. To fully identify a residue both the sequence ID and residue name are needed.
As of 2020, the PDB (organization) uses PDBx/mmCIF as the primary format. The PDB format is only a legacy one. The presentation of the insertion codes (and many other things) in the mmCIF format are rather confusing, so it's worth explaining.
The mmCIF format introduced double identifiers for chains, residues and atoms: author's ID and "label" ID. Both look the same and are easy to confuse. For example: two chains named B and A by the author are renamed as, respectively, A and B by the wwPDB; both names are used in parallel: one chain is B(A), the other A(B). Poor design choice at every step. The same with the sequence ID. Originally, we had two tags:

_atom_site.auth_seq_id -- author's sequence ID, possibly with an insertion code (100A),
_atom_site.label_seq_id -- numeric-only sequence ID, numbered with respect to the sequence of this molecule (SEQRES).

Then, as I heard, numeric insertion codes were spotted in the wild. (In a database big enough any mistake that is not automatically checked for will happen). This made a problem with translation from the PDB format to mmCIF. The sequence number 100 with IC 1 would be interpreted as sequence number 1001. The numeric ICs must have been remediated in the meantime, because they absent as of 2020. But at that time (10-20 years ago?), instead of correcting the entries, the decision was made to move the insertion code to a separate column (pdbx_PDB_ins_code).
The documentation was never updated, so the auth_seq_id is still documented to store the sequence ID, not just the sequence number.
Here is an example from 3B9F:
                             !!  I                                      nn
ATOM   770  C CD2 . LEU B 2  46  ? -23.416 21.041  -27.401 1.00 19.94 ? 60   LEU H CD2 1
ATOM   771  N N   . TYR B 2  47  A -21.191 16.717  -28.119 1.00 18.73 ? 60   TYR H N   1

The sequence IDs here are 60 and 60A. Since the meaning of the insertion code is not widely known, it landed next to label_seq_id which has nothing to do with it. So again: the insertion code is just a part of the author's sequence ID (originally, in auth_seq_id), and the letter A above (10th field) is simply an extension to the number 60 (17th field).


text:
Algorithmic investigation of viral proteins for vaccine design
I wish to investigate how a virus could be attenuated through targeted mutagenesis by a priori prediction of attenuating amino acid mutations.
In particular, I am seeking bioinformatic solutions where the protein sequence of an exisiting live vaccine strais has been compared against it wild-type counterparts to identify key attenuating mutations and the key attenuating mutation(s) can be deduced without wet-lab site directed mutagenesis.
I presume surface antigens, such as the spike protein in COVID-10 would be the proteins of interest, but I'm open to suggestions.
I am a computer scientist looking to investigate the mutational parameter space of viral protein sequence data. A key approach I am considering is to use a mutational correlation/covariance matrix and then remove convolution artefacts. I hope this will identify compensatory mutations that will maintain a local interaction for example within a fold, providing a starting point for identifying attenuating mutations.
This may have application for example to COVID-19 and the goal is to produce an efficacious vaccine.

COVID-19 is a complicated genome and the starting point is to get to grips with its fundamental biology. I strongly disagree with the statement bioinformatics has no role to play - quite the contrary you will not get anywhere without it right now.

There isn't a vaccine for any coronavirus, and your question is generally about targeted attentuation, which is a complex area.
The basic building blocks for any vaccine development is virological understanding of the proteins involved in pathogenesis. I will focus on covid-19 as an example here.
The majority of bioinformatics work is based around the spike protein as is wetlab vaccine development. This gene is immediately after the 20 kb of ORF1ab, including the RNA dependent RNA polymerase, its around 3.5 kb but it might be a bit bigger. The spike protein is a major B-cell (antibody) antigenic determinant and the cellular receptor sites are located within it. Your bioinformatics could ultimately focus on targeting these specific sites, in particular using antigenic (antibody) prediction algorithms.
First you need some basic biology of spike and hydrophobicity prediction is a good starting point such as TMHMM here and consider the SARS spike protein structure for example here. Once you've a handle on those structures antigenicity scanning would be a next phase analysis and possibly homology modelling. This is because COVID-19 and SARS are closely related, but there are key differences. You finally consider the full range of amino acid mutations in current epidemic via NCBI.
At this point you would target a specific vaccine design strategy and have a clearer understanding on the approach you wish to take. Good luck!


text:
using blastn to search [error]
so I'm fairly new to using linux and blast - have encountered an error in my work.
I have created a fasta database with the results in this directory path:
~/class/hw/db
I'm currently trying to run a search against my created database against query sequences in a file called seqs.fas.
My code is as follows:
blastn -db $db -query seqs.fas -out q1results.txt
However I keep getting this error:
Error: Too many positional arguments (1), the offending value: seqs.fas
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: seqs.fas
I'm not sure what I'm doing wrong here, and would appreciate any advice to fix this, thanks!

I remember having this problem in the past, it was encoding related. Rewriting the -s manually (the dash sign) should fix the issue, please see this.

Problem is a copy & past from the BLAST documentation, all dashes
  ('–') need to be replaced by hand at the command line.

Thanks for the comment - I've tried multiple ways to get rid of this, including typing out from hand, copying/pasting, and even copying the '-' that shows up when i type blastn -help, however I'm still getting the same error. Just to be clear, dashes ('-') refer to the 'minus' sign next to 0 right? (I'm a Mac user)

The only think that I can think of is the $db variable. Does it somehow contain whitespaces or any other character that the shell would not understand? If the tilde sign (~) is part of what $db corresponds to (like in your question), you can try changing it to the exact path of your database.

I've just tried ```blastn -db ~/5030/hw2/db -query seqs.fas -out q1results.txt
but that just straight up shows me a ```BLAST Database rror: no alias or index file found for nucleotide database in search path```

This is another error, I guess you point out to the wrong database path.

Just to be clear, when i'm doing -db, to which should I refer to: the directory where I have my database, or the entire DB path?


text:
Train a new basecalling model for nanopore?
Basecalling on nanopore requires a machine learning neural network model, for example here
A popular software from Nanopore for it is https://github.com/nanoporetech/taiyaki, which has many people forking it.
However, the project requires an existing neural model to get started. I don't seem to be able to find a project for training a brand new model for Nanopore.


text:
How to construct the phylogenetic tree of the whole genus?
I am confused about how to construct the phylogenetic tree of the streptomyces genus. The taxonomy of the streptomyces genus is mostly based on the MLST method, but now I need to use the gene data from NCBI to construct the phylogenetic tree of the streptomyces genus. Should I use the whole gene of the MLST's housekeeping gene to construct the tree? Or I should just use the loci to build the tree? Or the other method? I am a junior in university, have little knowledge of the filed., but my tutor asked me to do that. plus: Sorry for my awkward expression~

The MLST genes are not used to contruct a tree, because each pattern of alleles determines a predetermined genotype, e.g.  Staphlococcus aureus ST 22 is the MLST type concerning the MRSA epidemic. You could build a tree, but its only purpose is to visualise the ST, so its a bit pointless but okay at least you will not mistype something.
Performing a whole genome tree on Strep for the genus is complicated but doable and there are missing genes and stuff and if you just want typing then its not needed.

So if typing is your goal I would submit your MLST loci to the Steptomyces MLST database,
https://pubmlst.org/bigsdb?db=pubmlst_streptomyces_seqdef

Thanks!!! I've calculated the ANI of the streptomyces genus(only the type materials), and now I'm running the Orthomcl the get the core genes(the number of which is expected to exceed 1000), after that, I'll build the phylogenetic tree of the genus~

Could you kindly recommend some lessons or books about the phylogenetic filed for me? My major is horticulture and I've learned the basic lesson of bioinformatics and learned the Perl, Python, Linux, and R language. But I lack the professional knowledge of bioinformatics. Thanks for the help!!!!

At your level use Mega X its easy to use. There's quiet alot of info in phylogenetics theory throughout BioSE. If you get stuck you have to ask a separate question (Bio SE rules)

Thanks for your answer, I had an Internship at a bioinformatics company, so I've run the gene family project~ I'll study in the BioSE!!! But sorry, what is biose? I can't google it...

oh, I got that now, BioSE means biology searching engine...English isn't my mother tongue...

Sorry, Bioinformatics Stackexchange has quite a few threads on phylogenetics.


text:
Coronavirus RNA structures?
Is there anything known about the RNA structures of coronaviruses? More specifically - do they have any interesting known structures in the translatable region, like RRE of HIV or the double loops in flaviviruses?
Update
Here is a recent development on the side of structure prediction.

Excellent @vadim, thank you this is very informative. You can of course edit your question and then place your answer below, because you have essentially answered your own question. Also thank you for your question. When I've a minute I'll answer your 'RNA structure and vaccines' in more detail

Thank you for this suggestion. I am thinking to do a bit more reading, perhaps I will have a more complete answer.

To add a more complete answer: the current coronavirus is closely related to the SARS virus that caused the outbreak in 2004, and on which much research has been done.
Here is a general review of the coronavirus epidemiology, life cycle etc.
I haven't found yet any materials about the RNA structures in the translatable region, however the structures in the 3' and 5' terminal regions are well studied and are known to be crucial for initiating viral replication and for its packaging.
This recent article addresses in details the RNA structures in 5' terminal of human coronaviruses, whereas this older article describes the structures both in 3' and 5' terminals.
Finally, the interesting algorithmic developments for the link that I added to the question are described here.

huge thanks. This is a very interesting question and your answer is particulalrly good. The only point I'd like to add is compensatory mutations which are well known in viral RNA secondary structure.


text:
How to seperate the genus into different group based on th ANI calculation?
I've calculated the ANI of the streptomyces genus, but I'm confused about how to separate them into a different group, could anybody kindly tell me how to do that? Thanks!!!

Oh thats easy you use clustering analysis

ANI = average nucleotide identity.
This will give you a 2x2 matrix. You reconstruct a phylogeny using clustering, the most important being neighbor-joining, UPGMA is very occassionally used.
The statistic is clustering algorithm and its done by pairwise reconstruction of sister groups. Once a sister-group is defined (usually on nearest-neighor both in neighbor-joining and UPGMA) the genetic divergence from the hypothetical node is established either using the mean (UPGMA) or additive criteria (neighbor-joining). Okay enough theory ...
You get Mega X and find out the import format for a 2x2 distance matric. Mega X however will calculate this for you. You then go to the build tree menu heading and click "neighbor-joining" (or "distance"). The tree will very quickly automatically appear.

Thanks~ Wish you have a good day (or night???)  : )

But isn't the result of ANI is an NxN matrix? How to translate it into a 2x2 matrix? I used the dcast() funtion in R to transform the data, from A to B.
A:
Streptomyces_zinciresistens.fna Streptomyces_aidingensis.fna    78.6608 909     2493
Streptomyces_zinciresistens.fna Streptomyces_gilvigriseus.fna   78.4136 719     2493
...
B:
Kitasatospora.papulosa Streptomyces.achromogenes Streptomyces.adustus
Kitasatospora papulosa                  100.0000                   81.0354              81.1125

Okay, first up I don't use R, I use Python instead, the only thing I'll use R for is ggplot2. I've never worked on Strep. What we should do is arrange a time to chat on the phylogenetics debugging forum

Sorry, but where is the phylogenetics debugging forum? ... I come from a web restricted area, have little knowledge about the foreign web... CRY

Thank you!!! You are so nice !!!

??? Where is the invitation? Should I check my Gmail message box? I can't see it in stack exchange's message box

Hi @Tangli83 really sorry I was late, if you are around the same time tomorrow please let me know

Yes, I'll be around at the same time tomorrow, it's ok.

Just a minute I'm in a meeting

okkk, it'll be fine, i'll memorize the words(or recite words)


text:
processing multiple fastq files with cutadapt
I have DNA sample from 5 pools, having 25 fastq files each. I am running cutadapt to remove the primers using this command
cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG A7_S7_R1_001.fastq -o sb_A7_S7_R1_001.fastq --discard-untrimmed
However, this is too time consuming.  I tried my hand at writing shell script
for file in /dir/*
do
  cmd cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

But it is not running .
Can someone help.
Thanks

I assume you don't really have cmd in there. Put and ls in just to see which files are seen.

Try this one
for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

Also you have omitted --discard-untrimmed from the original command. If you want that as well along with everything to run in parallel:
for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG --discard-untrimmed "$file" -o "`basename -s .fastq $file`_out.fastq" &
done


text:
What does the yellow mean in this image from Virus Pathogen Resource?
What does the yellow mean in this image of a SARS (coronavirus) protease from the Virus Pathogen resource database? It's from https://www.viprbrc.org/brc/

Sulfur atoms are shown in yellow.
The molecular viewer that you use is JSMol (JMol ported to the web).
Atoms are colored by element: grey C, blue N, red O and yellow S.
If you wonder how other atoms would be colored, see JMol's Default element colors, by periodic table.

Here's the direct link: https://www.viprbrc.org/brc/structureOperation.spg?accession=1P9S&decorator=corona&context=1583712525175

Another question I had was this: https://bioinformatics.stackexchange.com/questions/11612/are-there-any-sulfur-dioxide-molecules-in-this-structure Thanks

Looks like they're using the common CPK method of colouring, it's not unique for J(S)Mol.


text:
How many single copy gene are needed to construct the phylogenetic tree?
I ran the orthomcl program with genomes of the whole streptomyces genus(from ncbi) as input, but I only got 22 orthologous groups(streptomyces's genome size vary from 2M to 15M, which shows the great difference among them). And the analyst from the bioinformatics company told me that my result wasn't reasonable if I try to construct the phylogenetic tree with the 22 orthologous groups I got. He also told me, I'd better use almost 1000 orthologous groups to build the tree. My question is, how to know my orthologous groups are enough to construct a tree?
The program was run during my summer vacation(a long time ago...), but I still can't realize why "1000 orthologous groups are needed". And the analyst form company just tell me that it was the experience number. But when I read some article about comparative genomics, I see some papers are used just about 10 single-copy genes to build the tree, why???

Thanks, your research objectives are complicated and difficult to answer in a question and answer format, even for a phylogeneticist with experience in bacteria. The other issue is language, because potential typing errors can change the meaning of what you are doing.  Lets keep trying, I'm sure we'll get there.

Okay I'm not really sure so I've bountied your question. Please supply as much information as possible  about your 22 contigs, lets see if someone else could anwer this.

The one thing that strikes me is that OrthoMCL was written under David Roos supervision, which really means it was written with Toxoplasma gondii in mind, i.e. a eukaryotic parasite with notable levels of gene duplication. They tested it against Plasmodium falciparum, but that isn't surprising because both are members of the Apicomplexa. Personally I would suggest looking at a tool more orientated towards bacteria, albeit I'd like to see the output of OrthoMCL.
Prokka is a tool that might give better results, here, which is designed for bacteria, but personally I would extract the orthologes/homologues directly if the primary purpose is to contruct a phylogeny, i.e. I would parse annotated genes which are not "hypothetical".

Just wanting to add that I've had good experiences with using Orthofinder on bacteria, and it's more user-friendly than I'd have thought possible for something that uses a command line (I'm leaving this as a comment because I've never used OrthoMCL).

I would say, 120 markers is enough for phylogenetic analysis.
see https://gtdb.ecogenomic.org/


text:
What is a proper way for random subsampling of metagenomic data?
Let's say we have a metagenomic sample that is paired-end FASTQ files including 10,000,000 DNA reads collected using shotgun sequencing.
How would one make a random subsample of the mentioned metagenomic sample with for example 1,000,000 reads? I know about the seqtk function that could be used in such cases, but my question is more about whether it is meaningful to represent a metagenomic sample with a small subsample generated using this way or it has to satisfy certain criteria, e.g. in terms of the number of reads and diversity to be a valid representative of the metagenomic sample?
Knowing that some species might be very rare in the sample, would it ever be proper to do subsampling on metagenomic samples?

Why do you want to subsample in the first place?

To make more samples to be used for machine learning

Ah ... data augmentation, yep thats a common approach. I personally can't quite see what teh problem is, if you can perform ML surely you can use a randomisation module, e.g. Python's random or sample?

I am not NGS expert, I known ML. Essentially you are performing data augmentation which is essential for ML. What I question is whether a replicate of e6 from a population of e7 would really be sufficient for ML because you need sample sizes (replicates in this case) of >100000.
The answer is really simple you bootstrap the data (sample with replacement), to obtain your replicates. Numpy has a random module that will do this,
p.random.seed(345)
pop = np.random.randint(0,10000000 , size=1000000)

Use pop to select from fastaq. The other way in is to define each fastq sequence within pop then simply,
sample = np.random.choice(pop, size=1000000)

Or a combination of the two. If you can do ML you can do Python coding.

Dear Michael, thank you for the answer, but what I am asking is not how to do random subsampling but the idea whether it makes sense to do it for metagenomic data where it may remove key information of rare taxons. I don't think that subsampling should be done with replacement anyway since the subsamples would not be considered independent for training and testing.

Ok, in that case I don't understand the rationale. Sampling is used heavily in machine learning, it is also used in statistical analysis such as Monte Carlo and  null distributions. Generally it is used to ensure the observation is an outlier (except for ML)


text:
Are there any Sulfur Dioxide molecules in this structure?
Does the protein structure shown in the JSmol link have any Sulfur Dioxide molecules?

In the picture, the Sulfur element is yellow and Oxygen is red (see here)
https://www.viprbrc.org/brc/structureOperation.spg?accession=1P9S&decorator=corona&context=1583712525175
Link to mmcif file (I have the VS Code plugin installed):
https://files.rcsb.org/view/1P9S.cif

The answer, as provided in the comments by gilleain is No. The commenter wanted to know why I asked this. The reason is I heard sulfur dioxide as silent speech when thinking about antibodies that’s why I asked.
I searched the file https://files.rcsb.org/view/1P9S.cif for the SO2 molecule and didn't find it. I got confirmation in the comments of the question.
Would anyone be interested in an app that takes a .cif file of a molecule and allows you to search the file for any type of element or molecule in it? Would that help anyone with their research?

Sorry, I do not see the advantage. All three PDB databases allow ligand search by name, residue name or SMILES in the browser. The PDBe has a decent API (https://www.ebi.ac.uk/pdbe/pdbe-rest-api), which gives info of which ligands are present, without downloading the structure. Furthermore, most instances of elements beyond the 6 of life are crystallisation artefacts, so it is way safer to search for enzymes that are expected to have more exotic elements, e.g. tungsten in a W-enzyme as seen in PDB:2W7Z.


text:
Identifying proteins in a concatenated protein sequences
I have a 120-concatenated protein sequence in a fasta file. I would like to split by proteins.
How can I identify each protein?

Do you have a multi-fasta file (A FASTA file with multiple sequences, each with its own header) or a mangled FASTA file with one sequence that is in reality 120 individual sequences concatenated?

A FASTA file with only one sequence.

In that case, how do you determine where one protein ends and another begins? How was this file generated?

That is what I want to know. The file was generated using GTDBtk, a toolkit programme that identifies 120 genes marker to perform taxonomic classifications. I am using this programme for another purpose. I need the sequence of each protein, but I am afraid I am not going to be able to do it, right?

sounds really tough unless you know how it was joined in the first place

Rerun GTDBtk so that it does not output all protein sequences into a single fasta file. I am not sure what is the correct usage for it, but typical bioinformatics software should not do it by default. Alternatively, look for non-amino acids characters in the output as possible delimiters. Or perhaps these are subsequences of fixed length, rather than entire proteins, in which cases split the sequence into 120 equal length substrings.

The only way I can think is to break the protein apart on methionine residues, blast each amino acid locus and concatenate any locus that does not score a query protein position of '1'. The theory is simple every protein starts with a methionine residue, however there are plenty of M residues within the protein that don't represent the beginning of the gene.
Breaking the protein apart on M would be easy, but then automating and parsing the blast would be fiddly. Obvioulsy you will need to code it and if its a one off would it be worth it?


text:
How to get the intersection of peaks after peak calling using MACS2?
I have following 5 narrow peak files after peak calling.
K14_peaks.narrowPeak
K15_peaks.narrowPeak
K16_peaks.narrowPeak
K3_peaks.narrowPeak
K8_peaks.narrowPeak

I need to get the intersection of peaks into a new file (The peaks that represent all these files). I used following code but I got higher number of peaks not a lower number.
bedtools intersect -a K14_peaks.narrowPeak -b K15_peaks.narrowPeak K16_peaks.narrowPeak K3_peaks.narrowPeak K8_peaks.narrowPeak -wo > intersected.bed

Can someone help me?
Thank you.

bedtools intersect with multiple files after -b performs pairwise intersections between -a and each file in -b.
It does not intersect all files simultaneously.
See here, from the Bedtools documentation. The green section only contains intervals from -a that also intersect at least one interval in any -b file.

What I think you're looking for is bedtools multiinter (documentation here is unfortunately blank).
bedtools multiinter -i input [input...] will produce a BED file that covers every position listed in your input files, followed by how many and which of your input samples overlap that given position.
You can use this BED file and create something like an UpSet plot to visualize how many bases/peaks are shared between any combination of your samples.

Thank you very much. Do you think is it ok do intersection one by one? I did it one by one. First intersecting the first two and generating the intersected file and then intersected file is again intersected with the third one. I got less peaks in this manner. Is it ok to do that?

It's not wrong to do it this way, no. But it is more difficult, takes longer, and you lose a lot of information in the process. If you use multiinter, you can later select that regions supported by all samples with awk '{if ($4 == x) print}' where x is the number of samples (5 in your case). And you haven't gotten rid of the sites that are supported by 3 or 4 samples, say.

Oh. I'll try that way. Thank you very much.


text:
Why there is extra probeID in microarry Raw data (CEL) as compare to data in Series Matrix File & SOFT formatted family file?
#Series Matrix File(s)
>data_matrix<- getGEO(filename="GSE115167_series_matrix.txt.gz")
>dim(exprs(data_matrix))
[1] 49395    12

#SOFT formatted family file(s)
>data_soft<-getGEO(filename="GSE115167_family.soft.gz")
>probesets <- Table(GPLList(data_soft)[[1]])$ID
>data_soft1<- do.call("cbind", lapply(GSMList(data_soft), function(x) {
 tab <- Table(x)
 mymatch <- match(probesets, tab$ID_REF)
 return(tab$VALUE[mymatch])}))
>data_soft1 <- apply(data_soft1, 2, function(x) {
 as.numeric(as.character(x))})
>rownames(data_soft1) <- probesets
>colnames(data_soft1) <- names(GSMList(data_soft))
>dim(data_soft1)
[1] 49395    12

#Raw data (GSE115167_RAW) .CEL file
>data_raw <- ReadAffy(widget=TRUE)
>data_raw<-rma(data_raw)
>dim(exprs(data_raw))
[1] 49495    12

There are 100 extra probe ID in raw data.Those 100 extra probe ID are not found in the
GPL15207 platform.
 >tail(exprs(data_matrix), 2)
 #shows probe ID
 AFFX-TrpnX-5_at
 AFFX-TrpnX-M_at

>tail(exprs(data_raw),2)
#shows probe ID
ERCC-00172-01_at
ERCC-00176-01_at

ERCC-00172-01_at, ERCC-00176-01_at probe ID are not found in GPL15207 platform. I found there are 100 of them in raw data. Can anyone tell what are those probe IDs?

Cross-posted: https://www.biostars.org/p/426426/

they are ercc spikeins ? https://www.thermofisher.com/order/catalog/product/4456740


text:
Pausing index calculation
I have Pol_II data for Ctrl vs Knockout condition with two replicates each. I want to calculate the pausing index for the Ctrl vs KO condition.
Here I have given tag directory as input to Homer.Now what I see is as such,for each biological replicate it calculates the pausing ratio.How can I calculate as a whole like 2 Ctrl vs 2 KO replicate?
I'm giving a sample of the output.
Any help or suggestion would be really appreciated.
dput(head(CTRL_PTPEM_Pol_II_pausing))
structure(list(`Transcript/RepeatID (cmd=analyzeRepeats.pl rna hg38 -strand both -count pausing -condenseGenes -d Control_Pol_II_new_set/ Control_Pol_II_set2/ PTPN6_KO_Pol_II_PTPEN/ PTPN6_KO_Pol_II_set2/ -rpkm)` = c("NM_173803",
"NM_001301039", "NR_132775", "NM_022114", "NM_001011666", "NR_160728"
), chr = c("chr16", "chr1", "chr3", "chr1", "chr7", "chr19"),
    start = c(15395704, 150272234, 13618331, 3069153, 28686083,
    22527485), end = c(15400754, 150277284, 13623381, 3074203,
    28691133, 22532535), strand = c("+", "+", "+", "+", "+",
    "-"), Length = c(250, 250, 250, 250, 250, 250), Copies = c(2,
    1, 1, 1, 1, 1), `Annotation/Divergence` = c("MPV17L|M-LPH|MLPH1|MLPH2|MPV17L1|-|16p13.11|protein-coding",
    "C1orf54|-|-|1q21.2|protein-coding", "SNORA93|-|-|3p25.1|snoRNA",
    "PRDM16|CMD1LL|KMT8F|LVNC8|MEL1|PFM13|-|1p36.32|protein-coding",
    "CREB5|CRE-BPA|CREB-5|CREBPA|-|7p15.1|protein-coding", "LOC105376917|-|-|19p12|ncRNA"
    ), `Control_Pol_II_new_set/ Pausing Ratio` = c(0.982, 1.832,
    1.114, 1.247, 3.476, 0.632), `Control_Pol_II_new_set/ Promoter Reads (-50,200)` = c(0.125,
    0.748, 0.249, 0.873, 0.499, 0), `Control_Pol_II_new_set/ GeneBody Reads (200,5000)` = c(0.13,
    0.338, 0.208, 0.669, 0.032, 0.091), `Control_Pol_II_set2/ Pausing Ratio` = c(0.686,
    1.131, 1.054, 0.816, 2.205, 1.005), `Control_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.1,
    0.7, 0.1, 0.3, 0.3, 0.1), `Control_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.203,
    0.604, 0.089, 0.396, 0.068, 0.099), `PTPN6_KO_Pol_II_PTPEN/ Pausing Ratio` = c(1.08,
    1.362, 1.108, 1.224, 1.943, 1.2), `PTPN6_KO_Pol_II_PTPEN/ Promoter Reads (-50,200)` = c(0.151,
    0.602, 0.151, 0.301, 0.452, 0.151), `PTPN6_KO_Pol_II_PTPEN/ GeneBody Reads (200,5000)` = c(0.125,
    0.392, 0.118, 0.212, 0.141, 0.094), `PTPN6_KO_Pol_II_set2/ Pausing Ratio` = c(1.652,
    0.909, 1.551, 1.395, 1.774, 0.584), `PTPN6_KO_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.429,
    0.536, 0.429, 0.536, 0.321, 0.107), `PTPN6_KO_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.206,
    0.603, 0.229, 0.346, 0.123, 0.279)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

but I would say the pausing index is something that is invented, commonly used, but with no clear indication of what is the expected distribution or how it can be subjected to test like you've mentioned.
I assume you would like to get your work published, so going by the papers with IF > 5 such as this, and this, and this, the clear strategy is use a boxplot to show the difference between conditions, and apply a statistical test such as mann whitney u or t.test to say, there is a difference.
Let's try this, sorry not very good in R,using your dataframe which i named x:
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

df = x %>% select(contains("pausing") | contains("Transcript")) %>%
rename_at(vars(1),substr,1,10) %>%
pivot_longer(-Transcript,names_to = "condition") %>%
mutate(condition=case_when(
grepl("Control",condition) ~ "Control",
grepl("PTPN6",condition) ~ "PTPN6_KO"
)) %>%
group_by(Transcript,condition) %>% summarize_all(mean)

In the data frame above, I have averaged the pausing index between replicates (not very good but then again, most papers have only 1 replicate), and group them by condition.
# A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   <chr>        <chr>     <dbl>
 1 NM_001011666 Control   2.84
 2 NM_001011666 PTPN6_KO  1.86
 3 NM_001301039 Control   1.48
 4 NM_001301039 PTPN6_KO  1.14
 5 NM_022114    Control   1.03
 6 NM_022114    PTPN6_KO  1.31
 7 NM_173803    Control   0.834
 8 NM_173803    PTPN6_KO  1.37
 9 NR_132775    Control   1.08
10 NR_132775    PTPN6_KO  1.33
11 NR_160728    Control   0.818
12 NR_160728    PTPN6_KO  0.892

Then we plot:
p <- ggboxplot(df, x = "condition", y = "value",
                color = "condition", palette =c("#c72e29","#016392","#be9c2e","#098154"),
                add = "jitter")
p+stat_compare_means(comparisons=list( c("Control", "PTPN6_KO")))

your response are life saver let me go through the links and will update you. One doubt can you just take average of the replicates thought of doing that but was not sure if that is right thing because homer might have some model to calculate it

you're welcome. well let me know if you find a way. if you just want to make a statement like it pauses more in condition a vs b, i think the above is definitely good enough

it's better to spend time on something more interesting like what's in the promoter or SUPER enhancer, than thinking about pausing index

"Warning message:
Computation failed in stat_signif():
not enough 'y' observations " i did get this as im using full data frame

damn i just copy pasted your code and fail to see that...............

so it basically calculates promoter/gene-body the pausing index I do see a significant difference but i read in the paper where they show like the difference in promoter vs gene-body. So can i comment from the above output if the pol-II is high at Promoter or gene-body PTPN6_KO compared to Control condition?

you can make a meta-gene plot, combined with this pausing index above, will tell you, most likely you have more pol-II at stuck at the promoter... i think it's pretty obvious, your phenotype


text:
Coronavirus phylogeny and evolution
This new article presents the phylogeny of the coronaviruses and the placement of the new coronavirus in it. This made me wonder of what we know about the coronavirus evolution, particularly where it comes to evading human immune defenses: it is not mutating as fast as HIV - indeed, the genomes responsible for outbreaks in various countries are almost identical. In the same time it is unlikely to undergo gene reshuffling typical of influenza, since coronaviruses have single RNA genome.

Very interesting, Prof Gorbalenya is part of the ICTV for Coronavirus so this paper makes sense.

FYI coronaviruses have proof reading activity, unlike HIV, its a huge genome, unlike HIV. Thus the point mutation rate is unlikely to be as rapid as HIV.
The phylogenies in this paper are not really about the evolution or epidemiology of the virus. The trees used for viral taxonomy and form part of the International Committee for the Taxonomy of Viruses (ICTV) assessment in officially naming the virus. The ICTV is the authority in naming a virus. The ICTV essentially called it SARS2 as opposed to the WHO who used the name COVID-19. Thats the point of this paper and will be massively cited.
The broader evolutionary epidemiology of "SARS2" is for further studies. The Lancet paper on 22nd February gave some indication of this (and I suspect you know about that paper).

Indeed, the paper I cited is not directly related to evolution. The nextstrain group that I also cited is normally working on the molecular evolution of viruses. However for now they are just tracking the data.


text:
RNAseq without depletion?
I have been asked by my supervisor to conduct an RNAseq without experimentally depleting the samples. However I am quite unsure regarding how much sequencing depth per sample I would require in order to map a decent percentage of the ORFs if I were to proceed this way.
Has anyone ever done RNAseq w/o depletion and with what results? Would you go on with it or should I make the case for depleting even though it's gonna be more expensive.
Thanks.

Hi and welcome, there's not enough information in your post to allow a comprehensive answer in my opinion. RNAseq of which species is quite key here.

Did they want you to do poly-A enrichment instead? What is the actual goal? What is the organism?

Hi, I just realized I did not mention I am working with plant endopythic bacteria, so no chances of using poly-A enrichment. The species is Acidovorax facilis, with the goal of assessing differential expression when exposed to different root exudates.


text:
REMC SNP annotations for transcriptome prediction - epiXcan
I have been trying to use a pipeline for predicting gene expression for a target gene (Although you can do it for many) - the name of this method is epiXcan. Here is a link to the paper:
https://www.nature.com/articles/s41467-019-11874-7#Sec22
One of the very first steps in the pipeline is obtaining SNP 'annotations' from an epigenome reference consortium known as the REMC. In the methods paper they say:
"We first prepare eQTL statistics (computed with MatrixEQTL4) and SNP annotations (extracted from REMC https://egg2.wustl.edu/roadmap/web_portal/). For each eQTL tissue, we use the matched REMC tissue to
extract the corresponding annotations".
I have completed the first step of obtaining eQTL statistics. However obtaining SNP annotations is where I am lost.
In my case, I am only trying to predict the expression of one gene (DRD4) and so my region of interest is the pre-frontal cortex. If you visit the website they provide though, it doesn't seem as though annotating a large list of SNPs - at least in the way they are presenting - is possible.
Can anyone provide some insight or direction on how to proceed? How does one obtain SNP annotations from the REMC? I have tried emailing the authors of the paper and they have not been responsive.
Thank you,

After much research I found out how to obtain the annotations.
You must go to the REMC website: https://egg2.wustl.edu/roadmap/web_portal/
Go to the chromatin-state-learning tab, and choose the n-chromatin state model desired - the core has 15 states, then 18, then 25 (imputed states).
Following this, you must go to the download URL link and, you will find a lot of files, primarily with the following coding information: E0*.. This corresponds to a particular tissue/cell type and so you must choose the one relevant to your research question. The way they are correspond to each tissue/cell type can be found in the following paper:
https://www.nature.com/articles/s41467-019-11874-7#Sec22


text:
Tools or softwares for clustering evaluation
I am using a few clustering algorithms as well as my own tool for protein sequences. Also, I have benchmarked clustering results (i.e. ground truth) to compare with the results. Is there any tool/software/script to generate evaluation metrics for the clustering results.
The clustering results/benchmarked are in the following format
1  Name1
1  Name2
2  Name3
2  Name4
2  Name5
.....

UPDATE: I am evaluating the clustering performance using scikit

from your other comment you have already found scikit, if your number of clusters/ground truth labels is reasonably low (say <50), take a look at the contingency matrix as well, it can give a very clear view of what is going on in your clustering https://scikit-learn.org/stable/modules/clustering.html#contingency-matrix

Thanks, @Pallie !! Yes, scikit looks good as some clustering papers used the same measures as well. The number of labels for my work is reasonably high (thousands and millions).

The definitive answer is no, that is why they are not so popular. There are metrics that can be used:

you can examine "minimum distance" as a criterion under the assumption that the lowest distance between taxa represents the largest homology (e.g. if it is sequence data), which would be translate to other biological scenarios.
resampling, e.g. bootstrapping, jack-knifing and examine the tree structure of each node against the biological plausibility of the answer being correct.

For both points a biological a priori, i.e. biological plausibility of assessing the tree structure is required to confirm/reject the clustering method.
Clustering algorithms are not like the probabilistic methods, where the highest likelihood is the singular criteria and of course they do give different answers. Its a brief explanation I know but hope it makes sense.

Thanks for the explanation. I am now trying  V-measure (https://www.aclweb.org/anthology/D07-1043.pdf) using sklearn (https://scikit-learn.org/stable/modules/generated/sklearn.metrics.v_measure_score.html)

I am evaluating the clustering performance using scikit

Okay, but thats a different question. You appear to be implying you have input output associations for machine learning, this information was not present in the question.

I have mentioned that I have both output and ground truth. I am sorry for the confusion.

No problem, its just ML is a different approach for verification involving e.g. ROC and accurary indices. If you can do ML, then its really good.

Thanks, @Michael !! I think I have to use ROC for further analysis. There are some drawbacks using scikit for my problem


text:
Prediction model for spread of coronavirus
A previous thread on Bioinformatics StackExchange citing an article by the Economist which references a cartoon of predictive outcomes by the CDC (Economist article on coronavirus) suggests a model in epidemiology using the "basic reproductive rate", commonly known as R0. Although I'm no bioinformatician, I imagine that one could tweak a model based on similar infections like SARS and HDN1. Yet, I found no models predicting or giving confidence intervals for the spread of coronavirus. I wonder if the the data is still too sparse, e.g. we don't know if it's seasonal (mentioned in the above thread) and we cannot infer the virulence.
Would such a model be reliable? If not, how many months of data would we need to reliably use a model predicting the spread?

The original R0 calculation was developed in the early 1900s to model malaria transmission and is known as the Ross-MacDonald model so your idea is over 100 years old. The model was very famous for identifying the optimal strategy for malarial control by targeting the adult mosquito rather than the mosquito larvae. What it said was reducing the mosquito population and reducing malaria transmission were two different issues. So the bednet strategy is based in this idea and particularly spraying the inside of houses with insecticide.
R0 was famously used in the 1980s to predict the HIV pandemic when WHO thought it would remain localised in Africa.
The problem with R0 and Coronaviruses is we don't know the asymptomatic population who could also transmit the virus.
This population is potentially very large. Epidemiologists will attempt extrapolation from cruise ship quarrentine, eg the Diamond Princess where everyone is tested to understand this transmission parameter and project this onto the rest of the planet. Clearly this could have limitations or inaccuracies. The other issue I am baffled by is seasonal transmission and whether this will occur with COVID-19.
One thing R0 of COVID-19 will identify is the latent period transmission which occurs before the symptoms, i.e. the incubation period. Essentially, people with cough/mild symptoms before they self isolate and become really ill with a lung infection will transmit the virus. In this context R0 will be/was very influential for understanding when to move between quarrentine and containment phases of the pandemic.

Have other models improved on this 100 year old model? It is now around 40 days since the quarantine of the Diamond Princess and around 2 weeks since they all disembarked; when would we know the asymptomatic population on that ship?

The 100 year old model is for malaria. The fundamental idea of a (R0) compartmental model which is solved by integration is very much alive.

The social stuff on the Diamond Princess I don't know about, I just know it is an ideal population to help solve the asymptomatic/transmission issue.

@miguelmorin there are many studies on this actually, but as Michael explains, the R0 values are all estimates and vary wildly. Just to give you an idea, this study, which is being used to inform the planning of the UK government (and probably others) posits an R0 of 2.4, while this one (which has not been peer reviewed, and which makes some pretty enormous assumptions; I wouldn't trust it much, personally) suggests R0 could be as high as 26.5.

To illustrate just how different those values are, if we assume R0=2.4, then if I am infected, I will infect an average of 2.4 other people. If they then do the same, after 10 steps, I will have ultimately infected 2.4^10 = ~6340 people. If R0=26.5, then after 10 steps, I will have infected 26.5^10 = ~1.7e10^14. That's 17 followed by 13 0s, one hundred and seventy trillion people.


text:
How to construct the phylogenetic tree of the "High Recombination" genus streptomyces?
I need to construct the phylogenetic tree of the whole genus of streptomyces. And my goal is to find out how my target strep species distributed in the whole strep genus and obtain the distribution pattern. Strep is a big genus, has 658 species(65% of them have the genome info), and my target species is about 25. But the traditional 16S method seems not helping when it comes to strep is a genetic recombination genus. And still lack of MLST scheme to distinguish my target species. I want to get the phylogeny for the whole genus, but I lack related genetics knowledge.
There is something I did before try to get the tree.

I ran the orthomcl program for 400 strep genomes which I got from NCBI, and the result is: I got 22 orthologues
I ran the fastani program for 200 strep genomes(type material) and got the ani result.

What can I do to achieve my goal? Thanks!

Hi @Tangli83, I will construct an answer in tomorrow evening (Saturday morning in China) based on our conversation

Oh, thanks,~ I thought my post didn't explain clearly.

that's fine I know your questions and research aims

You need bootstrap support using a model-based tree building algorhithm, via maximum likelihood (a few people use Bayes). The file format is relaxed phylip format (please submit a separate question if you have difficulties here - its a bit tricky).
I use standard RAxML here, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.
A robust maximum likelihood tree is,
raxmlHPC -m GTRGAMMA -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,
raxmlHPC -b 459 -m GTRGAMMA -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

This tree will bootstrap for 500 replications, however to start I would use 100 replications.
Make a consensus tree of the bootstraps,
raxmlHPC -f b -m GTRGAMMA -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)
The cluster you have access to is fine for the calculation, it will take aound 24 hours to complete one bootstrap calculation for one data set and obviously you need to parallelise you calculation across the 22 contigs.
Viewing the tree, FigTree (for Mac OSX) is easy.
Root
Rooting can be complicated because I don't really know your bacteria.

The recombination issue is more complicated, but I would construct 22 trees from your contigs and assess them for congruence. Panmixia is a concerned, which means too much recombination


text:
What Cellbase resources are 0-indexed, and which are 1-indexed?
The UCSC Genome Browser Team is clear about that the genome browser GUI indexing which is 1-based, closed interval. Genomic sequence retrieval via the UCSC REST API is 1-based, closed interval. All UCSC tables are indexed at 0-based, half-open intervals.
However, I can't seem to find a list of how resources are indexed at Cellbase. I know genomic sequence queries are 1-based and closed-interval, but are all the other tables? I could come up with tests for each resource I use at Cellbase, but I'd need a known true value for each test, which is impossible or impractical for most resources I need.
It's such a fundamental thing to know that there's probably a description somewhere, but I can't seem to find it. Or maybe there is a method to easily know, but I don't know that method. I've looked at:

the Cellbase website,
the Cellbase wiki (down for maintenance as of 3/11/20), and
I've read both the original NAR paper and
the recent BMC Bioinformatics paper on the pycellbase library (which I'm super-excited to try).

There is no mention of how nucleotides are counted. I'm sure I'm missing something; any help would be appreciated.


text:
polishing assembled genome to QV50 value
I am trying to assemble some genomes from pacbio reads using HGAP3 followed by quiver polishing. However, even after multiple rounds of quiver polishing, I cannot attain quality value 50 (QV50) for a genome.
QV50 seems to be recommended and I am able to get upto QV30 only. I have done polishing for almost 15 times already. I have tried polishing with different parameter changes like subread length cutoff and blasr maximum divergence.
Can anyone suggest me on how I can achieve an assembled genome of QV50 value?
Thank you

15 times is too many times. If you were going to get quiver to get you good QVs it would have happened after 2-3 rounds probably.
Some questions:

Why are you using quiver? Are your reads from RSII chemistry? Arrow is recommended for Sequel reads.
What is your dataset coverage? If you do not have very high coverage (I think >30X is usually a minimum for a eukaryotic genome) it may be that you are just not going to get exactly the values you want because assembling low coverage data is hard.
Have you tried other polishing tools? They might be able to get you better value. I am not an expert but some people like to use Racon or Apollo.
While I'm not an expert, looking into how the QVs are calculated can be helpful. I looked at the quiver docs and it looks like if the reads map ambiguously then QV is low. So if you for example have lots of haplotigs, then mapping might be ambiguous for a lot of reads, and therefore this will depress your QV values. This is therefore a problem that can come from the starting assembly and not just the polishing. You may want to go back and take another look at your starting assembly.


text:
Removing common variants in the 1000 genomes database from .vcf
I have 15 .vcf files. I need to remove `common variants in the 1000 genomes database' appearing in at least 0.5% of the population
Do you know from where I may start?
Thank you so much

Use VEP, ANNOVAR or snpEff to annotate your VCF file (I'd recommend combining your VCFs into a single file if they're single sample VCFs or are all comprised of samples from the same experiment/cohort). You can annotate the VCF file with 1000g (among a ton of other annotation sources). Once done, you can use bcftools view to subset the VCF as required.
Or, you can subset the standalone 1000g VCF and pick regions with your frequency threshold, then you can use that regions file to subset your VCF file(s).
Remember to left-align and normalize (and split multi-allelic records in) your VCF files so all comparisons can be done accounting for chr, pos, ref and alt fields.

Think removing everything < 0.05 from 1000g2015aug_all column of ANNOVAR annotation file helps with this

It's quite rude to add and accept only your own answer after people have attempted to help you, especially when your answer is incomplete and will not help others that may have the same question. I think I will stay away from your future questions and not attempt to help you, as you seem to do this a lot (evident from your pattern on biostars, where you also ignored others' help and only accepted your own half-baked answers)

@RamRS you are not answering this OP's question only, but many other future users'. Keep answering, time and number of votes will point to the right answer.


text:
are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?
I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:

MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.
mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default
FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.
GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.

The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length <500, which leaves a total of 21723 genesets tested.
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.

I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.

That's a really weird (and interesting) result. Do you know if each package came with its own library of genesets ? Or are they all using the same ?

Without code and details on how exactly you ran the tools (and with which input) this is hard to answer.

@Atpoint Although that is a perfectly valid observation, the main question is not about what could be happening and more about whether anyone has ran into similar results when comparing different methods. I edit the main post adding further details.

@dc37 Yes, the list of genesets was the same for all the tests, they are KEGG, GO:BP, GO:MF and GO:CC.

Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.
Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.
This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)
The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): Toward a gold standard for benchmarking gene set enrichment analysis.
Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.

Thank you, for some reason I was expecting that all the methods would agree on which genesets are differentially enriched to a certain degree.

Yes, I know it is somehow an expectation of my wet lab colleagues.