text:
Problems in creating desired phylogenetic tree with ggtree
I am working on haplotype data and want to make a tree out of haplogroups using ggtree. I have following data in newick format.
(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));

I am using following R code to plot a tree from this data,
> library("ggtree")
> tree <- read.tree(text = "(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));")
> p<-ggtree(tree)
> p+geom_tiplab()

And the displayed tree looks like this

I have following problems in this tree,

It gives a low quality image
It does not show the haplotype taxa names clearly
I want to arrange haplogroups i.e nodes in the desired according to the tree given below,

i.e B haplogroups with C haplotypes and with C, E named haplotypes etc.
I have tried to solve these problems but failed to find any solution. Is there a solution?

for the tiplab issue, pls refer to the FAQ.

The above given problems were solved using manual tree construction through MEGA. MEGA allows manual construction/manipulation of trees. The newick format obtained after manually constructing the desired tree using MEGA is:
((((((((((((T_M184,T_L131,T_PF5633),T),((R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),R)),((Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),Q)),(((((((O1a1_B384,O1b1_F2320),O1b2_P49),O2IMS_JST0213554),O2a1_F51),O_O2a2_F525),O3_M122),O)),(((N3,N1a1_M46),N1a2_L666),N)),(((((L_M20,L1_M22),L1a_M27),L1b_M317),L1c_M357),L)),(K,K_M9)),((((((((((((((((((((((J_M304,J1a_PF7257),J1a2a1a2_P58),J1a2a2_PF7264),J1a3_Z1828),J2_M172),J2a_L25),J2a1_L26),J2a1_M319),J2a1_Z1846),J2a1_Z387),J2a1_Z467),J2a1_Z6046),J2a1_Z6063),J2a1_Z6065),J2a1_Z7671),J2a1_Z7700),J2a2_PF5008),J2b_M241),J2b1_M205),J2b2_M241),J2b2_Z2456),J)),(((((((((((I_M170,I1_M253),I2a1_Isles),I2a1_S21825),I2a1a_M26),I2a1b3_L621),I2a2a_M223),I2a2b_L38),I2c2_Y16419),I),((H1_M69,H1a_M82),H)),((((((((((((G_M201,G1_M342),G2a1_Z6552),G2a2a_PF3147),G2a2b2b_PF3359),G2a2b1_M406),G2a2b2a1a_U1),G2a2b2a1b_L497),G2a2b2a1c_CTS342),G2a_L30),G2b_M3155),G2b_M283),G))),((E1b1a_V38,E1b1b_M215),(E2_M75,E))),(((((((C_RPS4Y711,C1a_V20),C2_M217),C2b1a_F3985),C2C1_Z1338),C),B2b_M112),B));


text:
No MQ tags in VCF files
To call minority variants in my Mtb sequences I'm using a pipeline of
fastq ->
    BWA MEM for mapping ->
    picard MarkDuplicates ->
    GATK3 for indel realignment and base recalibration ->
    lofreq for variant calling (https://csb5.github.io/lofreq/)

However the resulting VCF files don't have an MQ tag for me to filter on. Where in this process would mapping quality information be generated or lost so that it doesn't end up in the VCF file?

In GATK output for vcf, i think there are metrics that make use of MQ, for example MappingQuality, MappingQualityZero, MappingQualityRankSumTest. you can check it here https://gatk.broadinstitute.org/hc/en-us/articles/360037425531--Tool-Documentation-Index#VariantAnnotations

As @ATpoint explains in his comment, the mapping quality (MQ) applies to read mappings and, hence, is contained in BAM files, but not VCFs. The MQ values carry important information that ideally the variant caller should use when determining the variant quality (QUAL in the VCF) of individual variant loci.
In the LoFreq documentation it says that

Choice of mapper: Bowtie and BWA-SW
BWA-SW assigns very low mapping qualities to mapped reads. LoFreq incorporates mapping quality by default into its variant calling model, which doesn’t help in this scenario and will reduce its sensitivity. Consider disabling the use of mapping quality for lofreq call with -N.

So LoFreq is already aware of the MQ values and takes them into consideration during variant calling.
Since you are using BWA-MEM (and not BWA-SW) you should probably not be concerned about the issue they mention and trust their quality filtering.
If you are sceptic about the quality of mappings, though, you could still filter them out directly after mapping, e.g. using samtools view -q 20 ...

source: Stack Exchange

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

sort: votes


text:
What is tRNAscan-SE bit score?
I have used the command-line version of tRNAscan-SE by Lowe et al. (1997) on the reference genome sequence of C. elegans to predict tRNAs. The programme outputs a table with tRNA locations, amino acid, anti codon, intro locations, tRNA type (functional or pseudogene) along with a score. The score in my dataset ranges from 20.0 to 107.4. After reading the article Lowe et al (2016) describing this tool, I still do not understand how this score comes about and what it means. I understand that:

the higher the score, the more likely the tRNA is assigned to the correct isotype (amino acid)
the lower the score, the more likely the tRNA is a pseudogene

My questions are:

How is the score calculated?
What does the score mean?

The authors of the tool state that is is 99% sensitive and detects 1 false tRNA in 15 GB of DNA. Does this mean that I do not need to filter my results? What about the confidence in the type of tRNA and whether it is functional or pseudogenised?
Lowe et al (1997) tRNAscan-SE:  a  program  for  improved  detection  of transfer  RNA  genes  in  genomic  sequence. Nuc. Acids 25
Lowe et al (2016) tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes. Nuc Acids Web Server issue: W54–W57.

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


text:
Simulating 3' end tag-based scRNA-seq reads
Are there any tools that will simulate 3' end tag-based single-cell RNA-seq reads? For example, Drop-seq, 10X Chromium, CEL-seq?
There are tools that simulate scRNA-seq gene count data (e.g. Splatter), but I can't find anything that will simulate reads or even counts at the transcript level.

Looks like a couple options have become available:
Minnow: article code
scReadSim: preprint code


text:
Statistical approach to link DNA methylation with toxic element exposure and health outcome
I would be thankful to you if you can help with the statistical approach for case-control study to link DNA methylation (epic array) with toxic element exposure (arsenic) and health outcome (Cardiovascular diseases).
I do have data from 125 Cardiovascular diseases and 125 non-Cardiovascular diseases individuals matched with age, sex, and smoking.I do have phenotypes files in a csv format having sex, age, smoking, disease, samples id, iAs (arsenic) values in log format and then a separate files for methylation data.

Notes

The sample size is small.
Maybe structural equation modelling could be used to find CpG sites associated with arsenic exposure and disease condition as I have very small size samples
R is a preferred language.

This forum is not dedicated to general help requests and your post may be closed. Ask a precise bioinformatics question if you are encountering difficulties in your analysis, there is a lot of tutorials available on the internet.

What computer language are you going to utilize in this project?  Are you planning on an ML/AI approach?  That's a major question to ask yourself.  If you are planning on a statistical approach, I highly recommend R and then Python.  If your skills are advanced, Id try an AI/ML approach!  What databases/sources did you use to get the phenotype files in .csv format?  What data types are the csv values in?  Are they integers, decimals/floats, strings?  Should you consider a .tsv files format instead?  What's the purpose to your research question to arsenic and DNA methylation?

The question is part of a previous question here https://bioinformatics.stackexchange.com/questions/21823/mendelian-randomization . If they are taken together, plus comments it's answerable. I will make the appropriate edits shortly.

Having given some thought to this, I'd still use machine learning. I would attempt to augment the negative controls however. Is there a suitable published data set? R program is caret but it high barrier to entry without appropriate support.
250 samples is just about enough. If you can augment that to 350 with negatives ... should be enough. What you do is dump all your data ("features") into a single record, including the methylation data and use toxic element exposure as the training target. The output you want is called "feature selection", what you are wanting is methylation sites to carry more "weight" than any other feature. Then you start other analytics to see what it's associating with. Its not trivial getting it working and a sceptic would complain about the lack of significance in the results, but its a good starting point.
GWAS might work, just considering methylation vs toxic element exposure. SEM, I've never used it so I don't know.


text:
How to convert a phylogeny to a dendrogram in R
I want to convert my phylogeny into a dendrogram so I can use it with dendextend in R to produce a tanglegram. I have made some progress but I keep encountering errors, see below:
library(ape)
library(dendextend)

Tree <- rtree(10, rooted=F)

Tree <- read.tree(file="clipboard", text=NULL)
# test <- as.dendrogram(Tree)
## Error in ape::as.hclust.phylo(object) : the tree is not ultrametric

Tree.ultra <- chronos(Tree)
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not rooted

Tree.ultra$root.edge <- 0
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not binary

This thread seems to be the closest to answering my issues but the workflow is different and Im not sure how to address this binary issue. I have tried using as.hclust, as.hclust.phylo and hclust.

Please [edit] your question and give us a minimal reproducible example. Something that lets us test your code (and ours) and make sure our answers would work for you. A simple version of the tree you are using, for example, or even a way to generate random data that are similar to what you're working with.

Its not really clear to me (and maybe it is just me) what you mean? A dendrogram and a phylogeny are synonyms. Unless theres some suble difference I'm not aware of, so I'll let others enlighten me if thats the case.

Ive added a small example which incurs the same errors. @JoeHealey. Depends whom you ask about dendrogram vs phylogeny. Many say they are synonyms but classically a dendrogram is seen as a simple graphical representation of hierarchy vs phylogeny which is based on more sound and substantial, evidence and data.

If dendrogram is the same as cladogram than it contains just a tree topology, phylogeny is more general term, it frequently contains branch lengths corresponding to mutation rate * time. See this post.

I was googling and it seems that dendrogram is a term used for visualisation of hierarchical clustering, which could be analogous to NJ-tree, but then it is not just a topology. If you mean cladogram, can you use this word?

After searching high and low I have found an answer from this thread
Workflow goes like:
library(DECIPHER)
dend1 <- ReadDendrogram(file="clipboard")

# Alternatively to load from existing environment
Tree$node.label <- NULL # Need to remove node labels
dend1 <- ReadDendrogram(textConnection(write.tree(Tree)))

This object can then be used in the dendextend package and tanglegram function.

Does that give different results than the chronos function?

Contrary to what you say, you can read a character string from within R by passing a textConnection to ReadDendrogram.

Thanks @KonradRudolph I was not aware of that. Im only a novice to be honest. Ill change my answer above to include this.

Hi @KonradRudolph I have been trying the text connect as so: dend1 <- ReadDendrogram(textConnection(Tree)) but I get an error. Error in textConnection(Tree) : invalid 'text' argument

@AudileF Did you read the documentation? What’s Tree? A representation of the dendrogram in text format?

@AudileF not sure if you've figured it out already but what you should have done is dend1 <- ReadDendrogram(textConnection(write.tree(Tree)))

I hadn't actually @august thanks a mill I'll include it above now.

For your example to work you should replace the rtree function with the rcoal function, see example below.
library(ape)
library(dendextend)
Tree <- rtree(10, rooted=F)
is.ultrametric(Tree)
[1] FALSE
is.binary.tree(Tree)
[1] TRUE

Tree_rcoal <- rcoal(10, rooted=F)
is.ultrametric(Tree_rcoal)
[1] TRUE
is.binary.tree(Tree_rcoal)
[1] TRUE

Tree_rcoal2 <- rcoal(10, rooted=F)

tanglegram(Tree_rcoal, Tree_rcoal2)

To make use of a file in newick format, you can try to convert with the chronos function.
s <- "owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
is.ultrametric(tree.owls)
[1] FALSE
dendrogram <- chronos(tree.owls)

Setting initial dates...
Fitting in progress... get a first set of estimates
         Penalised log-lik = -24.04992
Optimising rates... dates... -24.04992
Optimising rates... dates... -24.04908

Done.
is.ultrametric(dendrogram)
[1] TRUE

Hi b.nota. Thats all well and good but I have an actual newick file I need convert to a dendrogram. rtree was to give an example for troubleshooting. Thanks though.

Yes, I understand. But you didn't mention any newick file in your original post. Check [link] https://stackoverflow.com/questions/7445684/how-to-convert-a-tree-to-a-dendrogram-in-r this post from years ago, I think the chronos function might be helpful in your case.

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


text:
How to isolate genes from whole genomes for phylogenetic tree analysis?
I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA [or] all ribosomal-protein-coding genes" https://www.biostars.org/p/1930/. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).

Hi Daniel, I would be cautious to build a phylogeny on 16s rRNA sequences as informative variability is essentially nil. I would suggest a good starting point being concatenated MLST genes, and work from there. You have the whole genome seems wrong to focus on the least variable gene.

Yes. Using 16S makes sense for distant species but not for closely related ones. it is extremely likely that there will be next to no variability at all.

If you haven't yet annotated the genomes, do so with Prokka as suggested in some of the other answers. Once you've got that, you can easily pull gene features out of a genbank by name using biopython, once you've identified a few sequences you want to use as the basis for your sequence typing.

There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:
for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses

Roary is my go to but you have no say in what come out in the alignment. If you want to select your gene content use blast as I suggest. Concatenate the output and use an alignment software.

@terdon I added a code snippet to do the analysis. It can be used to isolate genes (plural). If the OP just wants to build a phylogenetic tree, there's no need to isolate specific genes. With almost 500 genomes from a single species there's unlikely to me much resolution from analyses of a single gene or from a handful of MLST genes.

Extract desired gene sequences using standalone blast
Simply provide a reference database with your desired output.
Set up your command and away you go. You can set the search up with a for loop for a batch of sequences. Command may look like
for f in *.fasta; do
   f=$(basename $f .fasta)
   blastn \
   -outfmt "6 sseqid qseq %" \
   -query $f.fasta \
   -subject reference.fna \
   > out/$f.fas
done

Watch the output as blast will spit out the gene as detected in + or - sense. If you want to only gather positive sense use -strandoption.
The default output I have here is tab output which requires a few sed commands to make into fasta.
sed -i \
-e 's/\s*$//g' \
-e 's/^/>/g'  \
-e 's/\s\+/\n/g' \
*.fas

Online alignment servers are an easy way to align small datasets e.g. EBI

It might be worth mentioning that this approach work for species with no splicing (so it should be fine for the OP), but isn't very useful for eukaryotic organisms. Also, why would you run multiple blasts in a loop like this? Why not give a multifasta input file to blast and run it once?

In my experience the blast out uses the gene names from the reference. and does not append the original genome name. I go about a few QC/ parsing bits and then combine the sequences.

I'm not sure what my "desired gene sequences" are supposed to be. What is my query?

Your query is the genome sequences you have $f will insert the name for each in a loop. your desired gene sequences will be in your reference file. This could be the 16s or MLST gene sequences. Whatever you want to build a phylogeny of.

My question is more how can I find those MLST gene sequences in the first place?

A quick google brought this up http://bigsdb.pasteur.fr/klebsiella/klebsiella.html . Have a look about there is usually a file full of genes and alleles. In my experience one (allele) set of the 7 genes should suffice.

@DanielHarris take one MLST gene and use that as a query to find the rest.

@AudileF blast can take a multifasta file as a query and the results for each sequence in the input file will be shown in a separate Query= section. The names are indeed the names taken from the input file, but that shouldn't be a problem as long as they're all unique.

So use one of these MLST genes such as from this link http://bigsdb.pasteur.fr/perl/bigsdb/bigsdb.pl?db=pubmlst_klebsiella_seqdef_public&page=downloadAlleles as a query. Two questions now... Do I use the entire locus file, such as gapA has many >gapA_#'s or just one >gapA_1? Second, how does the output of blast then be used to build a phylogenetic tree? Blast output is far from multiple alignment format I think.

@terdon must admit I never know. Does this work in the case of contig assemblies?

@DanielHarris 1) use one copy of each of the 7 genes. 2) you are correct the output can be quite jumbled. I usually use the -strandoption to gather the + and - sense sequences. Split them by taxa name. Order sequences in each file by name. remove fasta headers and space. Create new fasta header of file name. Then concatenate them all. For simplicity you could use a online msa server like on ebi.  Long winded and confusing. Im not a master computational biologist terdon might have a better suggestion.

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


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

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

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

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


text:
Visually summarizing phylogenetic networks
After having worked with Bayesian phylogenetic tree inference for some time, I am now trying out a method to infer network phylogenies. To get an idea of what it infers, I'd like to find a way of summarizing the results.
The output of the inference is a NEXUS file containing a stream of networks encoded in eXtended Newick as trees with duplicate node names, like this:
(((((((Kiraman:0.028228,Kui:0.028228)S1:0.0797301)#H4:0.35229)#H2:23.4168)#H1:15.8027,(Kamang:29.3316,#H1:5.45455[&&NHX:gamma=0.36825381342449554])S11:10.3482)S12:1.27331,((Klon:1.16727,((#H4:0.524833[&&NHX:gamma=0.33697036713230577],Kafoa:0.632791)S6:0.0615575)#H3:0.472921[&&NHX:gamma=0.6642547809404831])S2:35.6082,(#H2:22.2107[&&NHX:gamma=0.7544913491697585],((A-Petleng:0.421975,(A-Takalelang:0.168615,(A-Ulaga:0.145472,A-Atimelang:0.145472)S4:0.0231422)S5:0.253361)S7:0.902144,#H3:0.629771)S3:21.3468)S8:14.1045)S9:4.17761)S10:54.5301):0;

￼
There is no formal restriction on the shape of the network, and hybridization edges have a weight 0<γ<1. How do visually or formally summarize a stack of several thousands of these networks, and in particular the backbone tree structure plus a glimpse of the hybridization events they display?
I have played with some tree summarizing approaches and tried to first infer a summary tree and then glue some cross-edges in there, but the way I came up with did not look reasonable.
I am now considering to calculate distances between leaves and using network visualization tools like SplitsTree or NeighborNet to get a visual representation of that, but how shoud I represent the uncertainty of the distances in even a single network (not to mention thousands of them) in the input for these algorithms?
Are there any good approaches to summarize phylogenetic networks of the shape given above? Or to consistently simplify them to a shape that can be easily summarized?

Firstly you treat this as two separate independent calculations.

The traditional approach is simply to produce a consensus phylogeny using your Bayesian method and assess the equivocal bootstraps/Bayesian probability.
Feed your data set (not the trees) into NeighborNet within SplitsTrees. There are statistics that classify the level of departure from a bifurcating phylogeny.

Examine both outputs for orthogonal support for incongruence.
Note Phylogenetic network statistics are a relatively soft area of phylogeny (they are not used as hard evidence).
Networks from Bayesian trees
The final approach (which is what you are asking about), is to feed your Bayesian trees into Dendroscope, this is a GUI that represents the uncertainty between different bifurcating phylogenies as networks and you can use it to understand the incongruence between different trees within your Bayesian MCMC, because the phylogenetic networks directly corresponds to the level of incongruence.
Dendroscope is a nice Java API and is very easy to use.


text:
How to measure euclidean distance between points with vtkDelaunay3D package?
I'm working with python vtkDelaunay3D (scipy.spatial.Delaunay) package for a special purpose. Normally the package is used for triangulation/tetrahedralisation.
How do I measure euclidean distance between points with this package?

Hi, could you provide some more details so that we can all better understand which sort of data you are working with?

I gave up using vtk package.Instead, I have decided to use scipy.spatial.delaunay which is a more robust package

@Sara If you can provide a small working example, maybe you could add your scipy-based solution as an answer.

@bli thanks but I got the answer and did that project one month ago

I think you'll find the answer here, https://stackoverflow.com/questions/20025784/how-to-visualize-3d-delaunay-triangulation-in-python


text:
A tool or webserver for building PSSM matrix
I have some protein sequences and I want to build a position-specific scoring matrix (PSSM) for them and then upload this PSSM to NCBI PSI-BLAST. I used CHAPS program for this pupose but uploading the output PSSM gave me an error in NCBI PSI-BLAST. Do you know any tool or webserver for getting PSSMs for a group of protein sequences which then can work in NCBI PSI-BLAST?
This is a part of my PSSM file:
PssmWithParameters ::= {
pssm {
    isProtein TRUE ,
    numRows 28 ,
    numColumns 131 ,
    byRow FALSE ,
    query
    seq {
        id {
            other {
                accession "WP_000208753" } } ,
        inst {
            repr raw ,
            mol aa ,
            length 131 ,
            seq-data
            ncbieaa "MTTKRKPYVRPMTSTWWKKLPFYRFYMLREGTAVPAVWFSIELIFGLFALKNGPEAW
AGFIDFLQNPVIVIINLITLAAALLHTKTWFELAPKAANIIVKDEKMGPEPIIKSLWAVTVVATIVILFVALYW" } } ,
    intermediateData {
        freqRatios {
            { 0, 10, 0 } ,
            { 564418841, 10, -10 } ,
            { 0, 10, 0 } ,
            { 11768571, 10, -9 } ,
            { 185838265, 10, -10 } ,
            { 31496547, 10, -9 } ,
            { 3872857, 10, -8 } ,
            { 291750464, 10, -10 } ,
            { 128450763, 10, -10 } ,
            { 759856221, 10, -10 } ,
            { 359173937, 10, -10 } ,
            { 179865517, 10, -9 } ,
            { 14537895, 10, -8 } ,
            { 212921456, 10, -10 } ,
            { 220554141, 10, -10 } ,
            { 368516324, 10, -10 } ,
            { 319343525, 10, -10 } ,
            { 426173953, 10, -10 } ,
            { 463659523, 10, -10 } ,
            { 817322186, 10, -10 } ,
            { 811693964, 10, -11 } ,
            { 0, 10, 0 } ,
            { 227810064, 10, -10 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 768325726, 10, -10 } ,
            { 0, 10, 0 } ,
            { 1425562, 10, -8 } ,
            { 372685285, 10, -10 } ,
            { 466727421, 10, -10 } ,
            { 185741084, 10, -10 } ,
            { 427328646, 10, -10 } ,
            { 122594965, 10, -10 } ,

Which error did it give you the NCBI PSI-BLAST? Could you paste here the image with the options used? What have you tried to overcome this problem?

The error was: Message ID#9 Error: PSSM processing error: Incomplete PSSM with Parameters provided on input: missing score matrix

So it seems that you need to provide the PSSM in certain format and CAPS didn't manage to create it in the expected format. Could you post the PSSM here (or in an accessible way) ? Now that I look in the help I couldn't find any reference to a specific format (just that it should be from their own output ...)

I do not exactly know how to post PSSM file here. Could you please guide me how to load a file here?

You can't but you can copy paste the content or add a link to dropbox or a similar server

@Llopis it should be fine posted here as code.

@Sara see formatting tools for help on formatting your posts. For things like this, you want to use code ({}) which uses a fixed-width font.

I think you could use standalone blast, one example here https://www.biostars.org/p/15482/

If you take a look at my answer in this BioStars post, you can generate a PSSM using AlignIO in Biopython:

I'm not 100% sure what you're trying to do with regards to your search in the first place, but I recommend using Sean Eddy's program Jackhammer it's similar to PSI-BLAST but uses Hammer HMM underpinnings.
As input it can take a single sequence, a multiple sequence alignment, or an existing HMM and then allow you to search in an iterative model building way.  You can also download an HMM too at any stage in the iterative search, as well as the results / alignment.

You can use biopython's code for generation of PSSM.
Here is the link for the source code which can help you.


text:
Find the date that sequencing was performed
Background I've obtained RNAseq datasets from ENCODE and I'm looking to identify external variables contributing to sequence variance among these samples.
External variables For example, clustering effects on the date the sequencing was performed.
Questions How do I find out the precise date which each sample went through the sequencing machine?

ENCODE databases are here on NCBI
The ENCODE website is here

What kind of data do you have? Do you have the count table or the original fastq files? Where did you get this files from (If you got it from someone could you ask him/her) ?

I downloaded them directly from the ENCODE website. I thought to ask here before e-mailing.

And why do you need the date? It is for batch finding? I don't know the metadata for the project but usually there documentation about how is this encoded

@Llopis: Yes, I'm looking to correct for variance associated with dates. Maybe I'm stretching it, but I'm really interested in parsing through the noise of the different datasets at all factor levels

Will you want to know which person did each part of the process? Or in which center, machine and lane it was done? I agree that this should be an objective, but I think it is unreasonable to ask this of this big projects, any way you can ask them. But I see there is a "Date released" filter value, so you can use that as a proxy

Yes, I agree-- for some parts it may be overkill. Date released corresponds to when the consortium released the data to the public. In fact, the sequencing lane and machines are in fact available in the library metadata! I thought maybe there is a way to get date information from a fastq file... I've read that dates are fingerprinted onto microarrays.


text:
How to find/build the evolutionary history of a protein from its sequence?
I'd like to build the evolutionary history of a protein, given its sequence. Namely, given a FASTA entry how can I build an evolutionary tree? Here is the 5wxy protein as an example:
>5WXY:A|PDBID|CHAIN|SEQUENCE
MGHHHHHHMMKTKLPILGVLGGMGPVVTAEFLKSIYEYNPFIDKEQESPNVIVFSFPSAPDRTGSIDSGKEREFIDFIQV
NLEHLNKLADCIVIGSCTAHYALPQIPENLKDKLISLIKIADQELQEYNEPTLLLASTGTYQKKLFQEGCTTADLIISLS
ESDQKLIHEMIYKVLKRGHDPLSILRDIEALLEKYNTRSYISGSTEFHLLTKSLKLKGIDSIKAIDPLSTIAQNFSQLII
KQAQVDLVTDCHQPSNPKSP

Do I need to do some queries to Entrez in order to build the evolution tree?
I saw this paper, but I am not sure how I can use it with a single mammal's protein primary sequence.

Question given a protein name, is it possible to look for its evolutionary history or build evolution tree using other databases beside PDB on top of bioinformatics algorithms.

You can’t build an evolutionary tree given a single data point (e.g. a single protein sequence). A tree relates multiple data points. Could you therefore please clarify what your input and desired output are?

@KonradRudolph, I would like to a protein from the PDB and find/generate its evolutionary tree.

Now you’ve just restated the question. But as I tried to explain, this doesn’t make sense.

@KonradRudolph I see what you are saying. Though how comes other people create such trees? for example (this paper) [http://www.pnas.org/content/108/20/8329.full].

The trees in that paper relate many different E. coli strains, which is exactly my point: rather than having a single data point, they are comparing several different ones.

@KonradRudolph so is there a database that associates a protein to its family etc.? I am not saying it's possible to get the whole tree from a single node. Of course it's not possible. My question is how can I do it using other resources given a PDB name.

Given a PDB identifier and sequence, one can use other databases like NCBI nucleotide database to find other related sequences (using psi-blast) to build a tree. All the resources you linked had a multiple comparison/alignment in its workflow, so you need other resources aside from PDB

@Llopis is there DB for amino acids database? Or an example for how to use what you have suggested? (Using NCBI/psi-blast)

There's uniprot, and you can use blast from protein to DNA. That is far away from your question. I'm just trying to help you to redirect your question for you to get answer. What have you tried/searched? I can spend some time helping people but if people don't spend time helping themselves I won't.

@Llopis what's blast? At the moment I try to understand what I need to know in order to find the evolutionary tree for a PDB entry.

Did you seach for it ? "blast bioinformatics" will lead you on the path

@Llopis, yes I ran it. What did you have in mind to do with it and http://www.uniprot.org?
The results: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=6KNTUB6D016

@0x90 What do you think was meant? You need to put a modicum of effort into things.

@DevonRyan I am asking for a general direction of what to do uniport offers various number of tools and I am not sure what should I do. I want to understand the idea of how to do it first rather to run bunch of tools without understanding.

@DevonRyan one of my problems is why it would be even possible to do evolution in reverse. I would expect evolution to be a irreversible process.

@0x90 If you are not sure what to do search for it, if even so you can't then solve your problem, then ask, but not the other way round, please. If you need more training you can take courses, there are some online courses very interesting

@Llopis I think what you guys gave me so far should give me some general idea to start with. If you have good idea for a course/reference on evolution in bioinformatics please let me know.

I am not informed of the online courses, but I'm sure you will find some, you can check back for opinions in the [chat] (if you want some).

The general procedure is:

Find the sequence for the same (or as similar as you can find (see mentions of blast in the comments above)) gene in other species. You can use uniprot or any other large sequence database for this.
You now have a large set of sequences related to yours, which you'll need to compare. Search the methods section of your favorite paper to see how they reconstructed likely evolutionary history from this, there are a number of packages out there.

My main concern is why these construction methods are plausible.

The exact same way we can infer what proto-indo-european words for things were...

lol how do linguists do this indeed? Also for languages you have ancient texts.

There you go, its Microcystis aeruginosa aspartate/glutamate racemase. Easy does it. If you want to know how you can do this in 15 seconds let me know.

Its one of the new features of NCBI's Blast

Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
Enter your sequence into the box (it doesn't accept PDB codes alone)
Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
Under the algorithm parameters enter "50" (default is too many)
Hit "Blast"
Once the search is complete at the top of the page are the hyperlinks:
"Other reports: Search Summary [Taxonomy reports] [Distance tree of results]"
Click on "Distance tree of results"
The following page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) .... here: https://www.ncbi.nlm.nih.gov/blast/treeview/treeView.cgi?request=page&blastRID=5A1ZBRU1014&queryID=lcl|Query_210080&entrezLim=&ex=&exl=&exh=&ns=50&screenWidth=1280&screenHeight=800
Click "Tool", "Download", "PDF" ...


text:
How to assess the quality of assembled .fasta genome files?
I have assembled 3 .fasta files from contigs infastq format of 3 different Homo sapiens.
I would like to see if the assembled genomes which are in .fasta format are in good quality or not. In addition, I would be happy to see heritages of genes from the maternal or the parental genomes.
More specifically, I have 3 .fasta files of a mother, father and daughter.

How can I see they are indeed the biological parents?
How can I study what genes were inherited from whom of the parents?
How can I assess the quality of the assembled genomes? Is there a standard way to do that?
Is there a simple frequency analysis of bases n-grams I can use?

Are these for humans? I imagine it'd be easier to work with VCF files than assembled fasta files.

@DevonRyan they are humans

What gene is this? Or it is the whole genome in a single fasta file? If the later case, if you explain the pipeline you use it might be easier to get an answer via VCF or tandem repeats.

@Llopis whole genome

To expand on the comment by Devon Ryan, what you really need are (A.) bam alignment files for each individual and (B.) VCF files for each individual showing all of their genomic variants (with respect to the reference human genome).  You mentioned that you had fastq files and that you converted them into fasta files.  Do you know if the bam files and VCF files are also available for download anywhere?
If the bam files are not available you will have generate them using sequence aligner software (for example: bwa-mem).  With the bam files you can calculate the average depth of coverage (and many other things).
If the VCF files are not available then you will have to generate them from the bam files using variant calling software (for example, the GATK package from the Broad Institute).
In terms of assessing the quality of the sequence runs, a commonly applied tool is name FastQC, but once again, it uses fastq formatted files as input, not fasta files.
The main use for fasta formatted files these days seems to be for compact storage of various reference genomes from different genome sequencing projects (using gzip or similar type of compression).
With the VCF files for each individual in hand you should be able to begin addressing your other questions regarding the frequency of shared variants, etc.

.bam and VCF aren't available. Is there any statistics that a valid genome should follow. Like frequencies of bases or ngram bases? Can you recommend of tools to generate VCF and .bam files from .fasta?


text:
changing blast parameters using NCBIWWW module
I have found a blog post with a script that I would like to use for my current research project: link
The script is incredibly fast and produces a smooth conservation plot. In the blog post, the author mentions that it would be totally possible to change the remote blast function to a local database. Instead of making a local database, I have opted to try to change the function to still use a remote blast, but with different parameters.
I found the following resource for learning about the NCBIWWW module: link
Using this, I rewrote the line in the script:
blast_handle = NCBIWWW.qblast(blast_method, "nr", search_gi)

to:
blast_handle = NCBIWWW.qblast(program="blastp",
                                      database="refseq_protein",
                                      sequence=search_gi,
                                      entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ",
                                      expect=20000.0,
                                      alignments=10000,
                                      descriptions=10000)

Strangely, I get the exact same graph as an output even with the different parameters. Am I doing something about this search incorrectly?
Thanks.

Update:
I got my BLAST to match the webpage. my current parameters are:
blast_handle = NCBIWWW.qblast(program="blastp", database="refseq_protein", sequence="P02686.3", word_size=6, expect=20000.0, hitlist_size=10000, gapcosts="11 1", matrix_name="BLOSUM62", filter="F", genetic_code=1, threshold=21, composition_based_statistics=2, alignments=10000, descriptions=10000, entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ")

With this new blast_handle, the script no longer works. I get the following error codes:

line 191, in
main(sys.argv[1])
line 27, in main
blast_rec = ncbi_manager.remote_blast(protein_gi, "blastp")
line 184, in remote_blast
return rec_it.next()
line 575, in parse
raise ValueError("Your XML file was empty")
ValueError: Your XML file was empty


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

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

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

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

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

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

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

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


text:
Is there a modern alignment tool tailored for transmembrane regions?
I am looking for a project or tool that allows programmatic pairwise alignments of proteins but that takes care with transmembrane regions of proteins. TM regions are traditionally too information poor for many tools to align accurately (Wong et al, 2010).

TM-aligner doesn't appear to have a program to download, and is only available via webserver. On the otherhand,
PRALINE has some optimisation for TM regions, however was released in 2008 and since then there has been a boom in TM knowledge from structures so does not benefit from that.Currently I am using PSI-BLAST. However,
PSI-Blast uses SEG to mask low complexity regions. This will most likely mean that TMHs are masked making the alignment of those regions fairly unreliable.

Primarily I am interesting in sequence alignments. But if anyone knows of a TM structural alignment, that could be good to know.
I think generally that counting idea works for an ideal helix, but a vast majority of TMHs break convention. I wouldn't expect to see a pattern based aligner that was that general. I'd be interested if anyone knows of a tool that juggles that though!
If I understand correctly, Matteo was saying that TMHs follow a simple patter (HHP/H,HHP/H...~7 times). This is not the case. There are too many examples of TMHs that do not match this pattern to assume this pattern would accurately align TMHs. As @M__ suggests, the number of TMP structures has not caught up to their soluble counterparts in the database. I've updated the question to explain why PSI-Blast might not be the best option.

Hmmmm I think you need to state an example, regarding the vast majority are not helices (suspect you are using GPI anchors as the definition). All viral TM proteins I know use helices.

Answer from @matteo-ferla converted from comment:
I assume you mean sequence (MSA), not structure alignment. Despite being a structural bioinformatician, I would slap a big caveat against structure based aligners. TCoffee:Espresso is probably the best structure-based aligner (browser and download and API), but the results are not as solid as some non-structural approaches such as training Psi-blast and aligning with that matrix. I am not familiar with MSA of TM regions, however, so this group may be an exception.
Given that helices are a repeat of 2 hydrophobic residues and 1.5 either polar (core-wards) or also hydrophobic (depending on how the helices pack), a pattern based aligner ought to do well from an evolutionary perspective. If you really want from a structure perspective and don't want to know about the cases where a loss of a turn at the top is compensated by the bottom, say, then TCoffee is definitely worth a try
Answer from @m converted from comment:
PSI-blast is a useful tool for deep alignment patterns, albeit it has its critics. I agree the vast majority of in vivo structures lack TM proteins, for different reasons


text:
What is a sensible forcefield choice for membrane proteins when using PDB2PQR?
I am generating PQR files for a membrane protein that is almost entirely buried in the membrane. The goal is to calculate the electrostatic charge across the surface of the protein. There are no lipids in the structure.
PDB2PQR has predefined options of amber, charmm, parse, tyl06, peoepb and swanson. It also allows user supplied forcefields.
The protein is an ion channel in the plasma membrane of a human. I'll be processing in memAPBS. However, I would greatly appreciate a "rule of thumb" answer as I plan to use the same method on TMPs from other organisms in various subcellular membranes.
What is a sensible forcefield choice for membrane proteins when generating PQR files for electrostatic analysis?


text:
How to allow "half windows" in biopython?
I am using Biopython to examine protein sequences including the Kyte and Doolittle hydrophobicity. However, for short sequences, this returns an empty value; the value is not allowed if the window does not have 20 values.
from Bio.SeqUtils.ProtParam import ProteinAnalysis

window_length=20
edge=1
kyte = {'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5, 'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, 'I': 4.5,
        'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6, 'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2}

full_sequence_analysis = ProteinAnalysis("CYGTGAVYLLGLGIGGFS")
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)

Is there are a built in Biopython option in this for "allow half windows/short windows" for incomplete windows? There is nothing I can see to suggest so in the docs.

@Chris_Rands This would fix the error. However, it would change the absolute values. Ideally, I would like to keep the window length at ±10 residues (up to 20 residues) if possible.

Have you tried odd values for window_length? https://web.expasy.org/protscale/ allows 3 5 7 ... 21

@nicolallias The issue would still remain. Even if I had 3 residue window length when I sample the final residue in the sequence, there are only 2 residues available for the window (the final residue, and the preceding residue), so no value is generated because 3 residues are needed for the windowed value.

You can choose between a window_size of 20 for sequences longer than 20 AAs, and the length of the sequence for those shorter::
seq = "CYGTGAVYLLGLGIGGFS"
window_length = min(len(seq), 20)

full_sequence_analysis = ProteinAnalysis(seq)
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)


text:
How to find which scales are available in SeqUtil of Biopython?
I am trying to analyse the hydrophobicity of a sequence using BioPython's "SeqUtil".
analysed_seq = ProteinAnalysis("AKVLGSLEALEPVHYEEKNWCE")
analysed_seq.protein_scale(Scale, WindowSize, Edge)

Known parameters

WindowSize refers to the total number of amino acids to be considered either side of the residue (an int variable).
Edge refers to the weighting of the scale further from the centre of the window (0-1). `

Unknown parameter

Scalerefers to which scale being used. I am interested in how to give theScale` variable a valid variable.

The Biopython docs simply point out that there are many scales to choose from and the source code doesn't seem to give any clues about scale options.
Where can I find a list of the scale are available and how to call them?

what do you mean? There is no scale argument to feed protein_scale() only def protein_scale(self, param_dict, window, edge=1.0): https://github.com/biopython/biopython/blob/master/Bio/SeqUtils/ProtParam.py#L211. I think param_dict is kd https://github.com/biopython/biopython/blob/master/Bio/SeqUtils/ProtParamData.py#L6

@Chris_Rands You've answered the question, thanks. I was looking for any predefined param_dict arguments and the GitHub source code spells them out clearly. It also looks like you can feed a scale straight in without having to edit ProtParamsData which satisfies me more.

Use this line of code:
print(analysed_seq.protein_scale(analysed_seq.get_amino_acids_percent(), window=9, edge=0.4))

Welcome to the forum.  Can you add some context as to what this line is doing and why it answers the question?


text:
Publicly available genome sequence database for viruses?
As a small introductory project, I want to compare genome sequences of  different strains of influenza virus.
What are the publicly available databases of influenza virus gene/genome sequences?

There area few different influenza virus database resources:

The Influenza Research Database (IRD) (a.k.a FluDB - based upon URL)

A NIAID Bioinformatics Resource Center or BRC which highly curates the data brought in and integrates it with numerous other relevant data types

The NCBI Influenza Virus Resource

A sub-project of the NCBI with data curated over and above the GenBank data that is part of the NCBI

The GISAID EpiFlu Database

A database of sequences from the Global Initiative on Sharing All Influenza Data. Has unique data from many countries but requires user agree to a data sharing policy.

The OpenFluDB

Former GISAID database that contains some sequence data that GenBank does not have.

For those who also may be interested in other virus databases, there are:

Virus Pathogen Resource (VIPR)

A companion portal to the IRD, which hosts curated and integrated data for most other NIAID A-C virus pathogens including (but not limited to) Ebola, Zika, Dengue, Enterovirus, and Hepatitis C

LANL HIV database

Los Alamos National Laboratory HIV database with HIV data and many useful tools for all virus bioinformatics

PaVE: Papilloma virus genome database (from quintik comment)

NIAID developed and maintained Papilloma virus bioinformatics portal

Disclaimer: I used to work for the IRD / VIPR and currently work for NIAID.

Glad that my upvote allowed you to expand your answer with more links :) I hear ViralZone http://viralzone.expasy.org/ is decent too (although I have not used it myself)

Thanks Chris! I agree Viral Zone is GREAT! I think it is one of the best virus education site around!

In addition to what others have suggested, I would also recommend PaVE as a resource. This is a curated database maintained by the NIAID and current holds over 300 papilloma virus genomes.

Hi quantik, thanks for your answer and welcome to Bioinformatics Stack Exchange. As a new user with reputation under 50 it's not possible yet for you to comment on other questions or answers (as @zx8754 has recommended), so in this case your answer was the most appropriate thing to do. If you want to improve on someone else's answer in the interim, feel free to edit your answer and incorporate a summary of the other answer(s) that you think are appropriate, making sure to include attribution where appropriate.

Thanks quantik...nice addition. A woeful oversight on my part as I am part of the group that develops PaVE. :-)

Influenza virus resource at NCBI or FluDB.

These are both included above. FluDB is actually the NIAID Influenza Research Database (IRD) Bioinformatics Resource Center (BRC) that has fludb.org as a URL.


text:
phangorn fasta file read error
Im trying to read a multifasta sequence file into phangorn I get the following error
tipseqfile <- ("abhi_seq/clean_dup_rem.fasta")
tipseq <- read.phyDat(tipseqfile,format="fasta")

my file

Error in read.phyDat(tipseqfile, format = "fasta") :
        sequences have different length

Now is it necessary to have same length of each fasta sequence? or Im doing something wrong

You have to perform an alignment first. Phangorn use phylogenetic estimation and this requires a precise alignment position to track homologous nucleotide mutation. Maximum likelihood theory will assess the probability of a given tree for each homologous nucleotide site, if they are not alignment it can't do it and will violate the assumptions of phylogeny.
If the sequences are different lengths, which it means both nucleotide and indel ('-') sites, it means there hasn't been or there is an error with the alignment. I use muscle these days having previously used clustalo. Its really easy to use and easily configures on a linux server and its super quick. If you get stuck let me know.
I had a look at your alignment and can verify it hasn't been aligned.

Once you've run your alignment, you may additionnally think about remove poorly conserved regions, e.g. with Gblocks. It might refine your phylogeny.

damn I was using the raw fasta file..got my error I have to use the alignment file


text:
Comparing phylogeny in R
So I want to compare the phylogeny created using two methods for example Maximum likelihood and maximum parsimony.Is there any way to compare the two phylogeny ?
I did read about phangorn but not sure if its the right R library for comparative analysis.
Any suggestion or helped would be highly appreciated
My data file
library(phangorn)
library(phytools)
library(dendextend)

data <- read.dna("abhi_seq/clean_dup_align_fast.fas", format = "fasta")
data
dat <- as.phyDat(data)

dm <- dist.ml(dat)
treeUPGMA <- upgma(dm)
treeNJ <- NJ(dm)

layout(matrix(c(1,2), 2, 1), height=c(1,2))
par(mar = c(0,0,2,0)+ 0.1)
plot(treeUPGMA, main="UPGMA")
plot(treeNJ, "phylogram", main="NJ")

dev.off()

parsimony(treeUPGMA, dat)
parsimony(treeNJ,dat)

tr.mp <- optim.parsimony(treeNJ, dat)

#tr.ml = optim.pml(treeNJ, dat)

fit <- pml(treeNJ, dat)
fit <- optim.pml(fit, rearrangement="NNI")

fit.ini <- pml(treeNJ, dat)
fit.ini

fit <- optim.pml(fit.ini, optNni=TRUE, optBf=TRUE, optQ=TRUE, optGamma=TRUE)
fit

tr.ml <- root(fit$tree,1)

tr.mp.ultra<-force.ultrametric(tr.mp)
tr.ml.ultra<-force.ultrametric(tr.ml)

is.ultrametric(tr.mp.ultra)
is.binary.tree(tr.mp.ultra)
is.rooted(tr.mp.ultra)

dd.ml.ultra<-as.dendrogram(tr.mp.ultra)

Error in ape::as.hclust.phylo(object) : the tree is not rooted

I ran into this error this error which says trees is not rooted
Comparing phylogeny code updated ** Working thanks to ***thomas duge de bernonville* thomas for putting codes together and fixing the errors
library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

a<-read.alignment("abhi_seq/clean_dup_align_fast.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a.rooted<-root(pars.a, outgroup="AAA64460", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))

mt <- modelTest(a.phydat, tree=upgma.a,multicore = TRUE,mc.cores=10)

#ml.a = pml(upgma.a,a.phydat)
#fitJC <- optim.pml(ml.a, TRUE)
#ml.a.rooted<-root(midpoint(fitJC$tree), outgroup="AAA64460", resolve.root=T)
#ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

#########################################################################

######################################################################33

#mt <- modelTest(dat, tree=tree, multicore=TRUE)
mt[order(mt$AICc),]# choose best model from the table according to AICc
bestmodel <- mt$Model[which.min(mt$AICc)]
env = attr(mt, "env")
fitStart = eval(get(bestmodel, env), env)
fit = optim.pml(fitStart, rearrangement = "stochastic",optGamma=TRUE, optInv=TRUE, model="GTR")#tree
bs=bootstrap.pml(fit, bs=25, optNni=TRUE, multicore=TRUE)

ml.a.rooted<-root(midpoint(fit$tree), outgroup="AAA64460", resolve.root=T)
ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

phangorn is a really powerful package for phylogenies. But to compare the trees, I think you may convert them into dendrograms and calculate a correlation measure, such as the Fowlkes-Mallows Index or a distance measure such as Baker’s Gamma Index. These can be easily computed using the dendextend R package (https://cran.r-project.org/web/packages/dendextend/vignettes/dendextend.html#correlation-measures).

"I think you may convert them into dendrograms and calculate a correlation measure" so you mean one for  Maximum likelihood and other for maximum parsimony, so can I do it using the Phangorn package? or I have to use something else

Yes, you can build the two trees and then convert them into dendrograms. Be careful, trees should be ultrametric for that.  You can convert them using the force.ultrametric function from phytools package followed by as.dendrogram from dendextend package.

I would let you know what i have done so far..So i have a multi fasta sequence which I did multiple sequence alignment using Muscle then i have an alignment file which i have used so far for distance based methods like NJ etc.The same alignmen file can be used in phytools to build Likelihood and parsimony tree ?

I think you can still use phangorn to construct the phylogenies.You just have to read your alignment with the read.alignment function of seqinr package, then convert it into a phyDat object, then calculate a distance with dist.ml. From that point, you may construct ML or parsimony trees (https://cran.r-project.org/web/packages/phangorn/vignettes/Trees.pdf). Once you've obtained the trees, convert them into dendrograms and calculate a correlation as described above.

data <- read.dna("clean_dup_align_fast.fas", format = "fasta")
data
dat <- as.phyDat(data)                                                                                d <- dist.alignment(dat, "identity") this is what i have done and I think as you mentioned Im doing it perhaps right need to see how to construct the ML or MP trees..

the ML and MP part I got in the vignette. How do I convert the trees into  dendrogram? this link https://bioinformatics.stackexchange.com/questions/923/how-to-convert-a-phylogeny-to-a-dendrogram-in-r  so this is what you are suggesting isn;t it?

from your phyDat object, you may estimate a distance eg d<-dist.ml(dat), then construct a rooted starting tree with UPGMA, tr.upgma<-upgma(d) then from this starting tree construct either ML or MP following the phangorn vignette. Once you get the bootstrapped trees, set them as ultrametric tr.ml.ultra<-force.ultrametric(tr.ml) and convert them as dendrogram dd.ml.ultra<-as.dendrogram(tr.ml.ultra). Please load the appropriate libraries to do that, and to have a look at the parameters that can be tuned at each step.

I have updated my question with code and the data file and the error i ran into when Im trying to covert it into dendrogram

I think you should use the UPGMA tree as a starting tree, both for the MP and ML, because I think it produces rooted trees, in contrast to NJ.

okay let me try that and update you i did try the UPGMA im still getting the same error "Error in ape::as.hclust.phylo(object) : the tree is not rooted
"

Let us continue this discussion in chat.

Another 2 tree compare tools: ggtree and blatic 3

The tree building is usually used IQ-Tree, MrBayes, RAxML, and beast to get a robust maximum likelihood or maximum-credibility tree. This process is time-consuming and computational resource consuming. If I want to use the dendextend to compare the phylogenetic trees which are built by the methods mentioned above rather than in R, can I do it with dendextend?

Bootstrapping I agree correlation is one way in, but the classic way to do this is via bootstrapping, i.e. resampling the alignment with replacement and making a consensus phylogeny of of 100 to 1000 replicates. Values exceeding 75 to 80% define node robustness and you can calculate incongruence from there.
Phangorn certainly does bootstrapping, but the efficiency of the algorhithm for maximum likelihood (ML) I don't know (might take a very long time), just think how long one tree takes and times it by 100 (or even 1000).
If you just performed a correlation measure there is a risk the differences are not robust, but if you've got loads of taxa you could bootstrap collapse any value less than the robustness threshold into a polytomy and then perform a correlation method. Its a bit complicated, but it would work.
Bootstrapping parsimony isn't hard nor computationally expensive, bootstrapping maximum likelihood sometimes requires specific efficient algorhithms notably RAxML.
Non-R The other way in is just dump the trees into Dendroscope3, okay its not R but will work. Dendroscrope is hard for a non-tree person to read, however it will flag the differences for a specialist. Dendroscope will make a network tree for the topological differences, in this case between methods, i.e. non-bifurcating. You really then usually have to present the differences as a "mirror tree", i.e. parsimony and ML side by side, because usually non-tree people wonder what a Dendroscope output means. Technically there is nothing wrong with describing incongruence via networks.
Basically, the heart-beat is bootstrapping.
Generally speaking parsimony will give the same basic output as ML, except when you have some taxa undergoing rapid evolution against their sister group. This does happen in pathogen evolution and parsimony falls over into a phenomina known as 'long-branch attraction'.

The advantage of the correlation measure will be the metrics which provide a more objective evaluation. Michael's recommandation with bootstrapping is really good, especially using collapsed branches.

If you only have to compare two phylogenies (as opposed to more), I would suggest that an alternative to using R would be using iqtree . One of its parameters (-rf) can calculate the Robinson-Foulds metric between two trees.
(you can also use iqtree to calculate the phylogenies in the first place)
It is command-line based, but can be easily installed with conda. Having anaconda installed in general makes bioinformatics work easier, I feel.
Alternatively, differences in trees can be visualized (without numbers) in Dendroscope, or Splitstree. They have different strategies for it. With Dendroscope it is easier to compare side-by-side and, if having precisely the same taxon names, to draw lines between taxa in different trees. Splitstree can accept an input of multiple trees (File>Tools>Load Multiple trees) and build a network (which can be a way to visualize the consensus of the trees).

Hi @Laura, SplitsTree is very cool, but by default will not work here because classically it will identify incongruence within the sequence itself (not just between trees). There may be options to disable this (if so please post them), but by default thats what it will attempt.

@Michael Thank you for pointing that out. I meant Splitstree starting from trees, not from sequence (File>Tools>Load Multiple trees). I've used it before for visualizing bootstraping results, but technically it should work on fewer trees too?


text:
Annotation label not matching to the color assigned in complex heatmap
The issue im facing is I'm assignig the colour to the annotation label but when it is rendered the label in my metadata does't match with the final output not sure what exactly is the problem.
The code
data file  ka_ks and metadata
As an example this sample Q6Q1S2 in my metadata is Alpha but in the figure it is labelled under Gamma
Any suggestion or help would be really appreciated
    data <- read.csv("ka_ks.txt",header = TRUE,row.names = 1,sep = ",")
head(data)

#metadata <- read.csv("NUC_META_AMINO_heatmap_metadata_sans_nucleotide.txt",sep = ",")
metadata <- read.csv('vir_meta.txt',sep = ",")
#head(data)
head(metadata)
#names(metadata)[3] = "Host"

#Set annotation
ann <- metadata[,c(2,3)]
colnames(ann) <- c("Class","Host")
#rownames(ann) = metadata$Sample
ann
#ann = metadata
colours <- list("Host"=c("Human SARS coronavirus 2"="green4",
                         "Human SARS coronavirus"="darkgreen",
                         "Bat coronavirus"="wheat4",
                         "Murine coronavirus"="turquoise","MERS"="orange",
                         "Bovine coronavirus"="coral",
                         "Rat coronavirus"="saddle brown",
                         "Human coronavirus"="dark green",
                         "Porcine hemagglutinating encephalomyelitis virus"="maroon",
                         "Avian infectious bronchitis virus"="dark slate gray",
                         "Porcine respiratory coronavirus"="slategrey",
                         "Canine coronavirus"="purple",
                         "Porcine transmissible gastroenteritis coronavirus"="pink",
                         "Feline coronavirus"="gold",
                         "Porcine epidemic diarrhea virus"="brown"),
                "Class"= c("Alpha" = "blue" ,"Beta" ="black","Gamma" ="red"))

heat <- scale((data))
head(heat)

colAnn <- HeatmapAnnotation(df=ann, which="col",col=colours, annotation_width=unit(c(1, 4), "cm"), gap=unit(1, "mm"))
myCol1 <- colorRampPalette(c("navyblue", "white", "red"))(100)
myBreaks1 <- seq(-4, 4, length.out=100)

#Add column boxplot
boxplotCol <- HeatmapAnnotation(
  boxplot = anno_boxplot(data.matrix(heat),
                         border = TRUE,
                         gp = gpar(fill = '#CCCCCC'),
                         pch = '.',
                         size = unit(2, 'mm'),
                         axis = TRUE,
                         axis_side = 'left',
                         axis_gp = gpar(fontsize = 12)),
  annotation_width = unit(c(1, 5.0), 'cm'),
  which = 'col')

hmap1 <- Heatmap(heat,
                 name = 'ks_score',
                 #name="Z-score",
                 col=colorRamp2(myBreaks1, myCol1),
                 heatmap_legend_param=list(color_bar="continuous",
                                           legend_direction="vertical", legend_width=unit(5,"cm"),
                                           title_position="lefttop", title_gp=gpar(fontsize=10, fontface="bold")),

                 #Split heatmap rows by gene family
                 #split=df$Family,
                 #split = NULL,

                 #Row annotation configurations
                 cluster_rows=TRUE,
                 show_row_dend=TRUE,
                 row_dend_side = "left",
                 #row_title="Transcript", #overridden by 'split' it seems
                 row_title_side="left",
                 row_title_gp=gpar(fontsize=30,fontface="bold"),
                 row_names_gp = gpar(fontsize = 10),
                 show_row_names=TRUE,
                 row_names_side="right",
                 row_title_rot=0,

                 column_dend_reorder = TRUE,
                 #row_dend_reorder = TRUE,

                 #Column annotation configuratiions
                 cluster_columns=TRUE,
                 show_column_dend=FALSE,
                 column_title="Samples",
                 column_title_side="top",
                 column_title_gp=gpar(fontsize=20, fontface="bold"),
                 column_title_rot=0,
                 column_names_gp = gpar(fontsize = 10, fontface="bold"),
                 show_column_names=TRUE,

                 #Dendrogram configurations: columns
                 #clustering_distance_columns="euclidean",
                 clustering_method_columns="ward.D2",
                 column_dend_height=unit(10,"mm"),

                 #Dendrogram configurations: rows
                 #clustering_distance_rows="euclidean",
                 clustering_method_rows="ward.D2",
                 row_dend_width=unit(15,"mm"),

                 #top_annotation_height=unit(0.5,"cm"),
                 #top_annotation=ColAnn,
                 top_annotation=colAnn)

#bottom_annotation_height=unit(2, "cm"),
#bottom_annotation=boxplotCol)

draw(hmap1 ,heatmap_legend_side ='right',annotation_legend_side = 'right')

fixed it ..the orders between my metadata and the data matrix was not matching


text:
Convert MeSH term to drugbank ID
I have set of annotated term of drugs extracted from PubTator those drugs are annotated in form of MeSH terms. I would like to map them into drugbank ID. I do search the the MeSH term i find the reference to the drug and link to the drugbank. But not any database as such. How do i map the MeSH term to the drugbank ID.?
My small dataset.
Terms         MESH
1        ara-C MESH:D003561
2 mitoxantrone MESH:D008942
3    etoposide MESH:D005047
4  fludarabine MESH:C024352
5   gemtuzumab MESH:C406061
6   ozogamicin            -

As an example of one of the drug I would give this drug "gemtuzumab"
DB00056 Gemtuzumab ozogamicin   Myeloid cell surface antigen CD33   target

So i want to map the above drugs into drugbank ID and add the information such as the class which is here target or it can be a transporter or enzyme etc.
The issue is there is no common ID between the drug MEsH and the drug-bank ID.
Any suggestion or help would be highly appreciated if anyone can cite me some resources.

PubTrator description from the website
PubTator Central (PTC) is a Web-based system providing automatic annotations of biomedical concepts such as genes and mutations in PubMed abstracts and PMC full-text articles.


text:
Is there a standard way to clean a PDB file and re-number its residues?
Is there a pre-existing tool which will tidy up the numbering of a PDB file?
Firstly, I would like to re-number the residues on inserts to make the icode an actual residue in the chain (by that I mean it's own number, shifting everything after it, so 45A becomes 46).
Secondly, I would also like to re-number so that each residue is one after another, as there are skips in the numbering, in spite of the two residues being physically next to each other.
I could do this in-house with a python script, but I thought I would ask specifically if there is a library or pre-existing tool?
Worst case scenario I have found this Perl script that does similar.
Thanks.

"make the icode an actual residue in the chain" -- something is confused here. Insertion code is part of the residue id, together with the sequence number.

@marcin Yes  sorry you're right it is, but for example when using the biopython PDB to FASTA parser it is not recognized and the chain is truncated from beyond there. So either I would need to ignore it (as I am at the moment) or change it from 1, 2, 3, 3A, 4 to 1, 2, 3, 4, 5.

@TW93 You almost certainly want to either ignore them, or deal with them separately.  They usually represent different wildtypes/isoforms.  E.g. the same protein within different species which contains some minor sequence alternations!  If you include the same residue twice or more, obviously, that will cause steric clashes and mess up any calculations that you run!

I think there is a confusion here between insertion code and alternative location ID. Insertion codes allow for unified numbering between proteins of slightly different sequences, but everything with a unique residue ID (residue number and insertion code) is a separate residue in the chain and cannot be ignored. The question asked by the OP is hence valid. You wouldn't want to have multiple residues with the same residue ID but different alternative location IDs renumbered, but that wasn't the question.

Several years ago we have written a tool for this: http://kodomo.fbb.msu.ru/~avkitex/projects/pdbparser.html
https://github.com/avkitex/pdbparserTools
The description is only in Russian but you can contact the author.

Could you explain how this program solves the question? Adding an example will be helpful to all users

The pdb-tools suite has pdb_delinsertion and pdb_reres for this task. Something like:
pdb_delinsertion 1CTF.pdb | pdb_reres > out.pdb

should work.


text:
What is the current state-of-the-art in assembling hybrid transcriptomes?
We are considering attempting de novo assembly of a species transcriptomes (i.e. without a reference genome) using the combined NGS outputs of Iso-seq and Illumina.
One example I saw (Li et al 2017), used the standard PacBio tools to assemble a transcriptome, followed by correcting the assembled sequence using a tool called proovread, followed by cd-hit-est and then cogent.
Would this be considered a reasonable pipeline or is there an alterative recommendation?

My recommendation would be to run Trinity using the --long_reads option, which allows you to provide an error-corrected fastq file for anchoring reads to a transcript isoform. The reads are then clustered per-isoform, and assembled, in a similar fashion to what's done with genome-guided Trinity:
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity
Trinity --seqType fq --max_memory 50G --long_reads corrected_pacbio_reads.fa --left reads_1.fq  --right reads_2.fq --CPU 6

Disclaimer: I have contributed code to Trinity at some time in the distant past


text:
Which sequence alignment tools support codon alignment?
Sometimes it useful to perform a nucleotide protein coding gene sequence alignment based on codons, not on individual nucleotides. For example for further codon model analysis it is important to have full codons.
A widely used approach here is to perform a protein sequence alignment first and then impose this alignment to the nucleotide sequences using PAL2NAL, CodonAlign or something similar.
This is how transAlign or GUIDANCE (in codon mode) work.
The problem here is that you are discarding part of the information which could be potentially used for the sequence alignment. E.g. if you have slowly evolving low-complexity region adjacent to a quickly evolving one, the amino acid induced alignment could be wrong, while incorporating nucleotide sequence potentially allows to make the alignment more accurate.
I'm aware of two programs which can do true codon alignment. First, PRANK has a dedicated codon model, but it is rather slow and using it is overkill for certain problems. Second, Sequence Manipulation Suite can perform codon alignments, but only for a pair of sequences; also it's javascript based, therefore it is hard to run it for a large number of sequences.
Can you recommend any software for multiple codon sequence alignment? Preferably available for offline use.

Why not translate the sequence and perform amino acid alignment? True, you “discard” the codon isoacceptor identity but that can be mapped back trivially if needed.

As I write in the question, this way you loose some information. E.g. you can produce an incorrect amino acid alignment for low-complexity region, which you otherwise would be able to reconstruct MSA properly.
In principle I agree that this doesn't happen very often.

There's blastx and tblastx but I don't think they do true codon alignment.

@IakovDavydov how does a codon alignment 'loose some information' over a simple DNA-DNA alignment? It's easier to find homologous residues using proteins (instead of DNA), and using that information to align DNA triplets can only be beneficial.

@PejoPhylo I was not claiming that codon alignment looses some information. I was saying that when converting codons to amino acid some information is lost.

@IakovDavydov I'm sorry, but then I don't understand what you originally asked in this question. Couldn't you use PAL2NAL for your (multiple) codon alignment? It supports offline use

@PejoPhylo I don't think I can make it much clearer here than in the original question. But briefly: with PAL2NAL you loose some information, which might be important for some alignments. Most of the time that's not important, but if your dN is comparable with dS and GC-content stays constant, that could make a difference. This paper provides some details: https://doi.org/10.1093/molbev/msr272

@IakovDavydov thanks, I'll give it a read.

@IakovDavydov I don't think that your original question is correct.  If one is interested in the functional aspect (the protein), then any codon based aligner will work, for example MEGA.  By definition, the amino acids are being aligned. The question about the low complexity, slowly evolving region is mute. If you are interested in the codon differences and drift, or evolution of the codons, then yes - but, I think that using nt aligned by codon is done more from the perspective of bundling the nt into triplets.

I don't know of any transcript-to-transcript aligners that are able to do this, but LAST can align transcript queries to protein reference sequences using a specified frameshift cost. Here's the specific documentation for that option:

-F COST
Align DNA queries to protein reference sequences, using the specified
  frameshift cost. A value of 15 seems to be reasonable. (As a special
  case, -F0 means DNA-versus-protein alignment without frameshifts,
  which is faster.) The output looks like this:
a score=108 s prot 2  40 + 649
FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD s dna  8 117 + 999
FFLQ-IKLWDP\STPH*IVSSP/PSDLISAHTLCPRMKSQDN

The \ indicates a forward shift by one nucleotide, and the / indicates
  a reverse shift by one nucleotide. The * indicates a stop codon. The
  same alignment in tabular format looks like this:
108 prot 2 40 + 649 dna 8 117 + 999 4,1:0,6,0:1,10,0:-1,19

The "-1" indicates the reverse frameshift.

I sent an email to the LAST mailing list about adding a frameshift penalty for transcript-to-transcript matching; I've been pleasantly surprised with the requested features that Martin Frith has added to LAST in the past. Unfortunately, in this case the problem is too difficult to sort out due to all the possible combinations that could happen, so it's unlikely to be implemented in LAST in the forseeable future (unless someone else writes that code).

Try MACSE v2 (https://academic.oup.com/mbe/article/35/10/2582/5079334) will align multiple protein-coding nucleotide sequences based on their amino acid translation while allowing for the occurrence of frameshifts

Thanks @user90!

I had no idea! thanks, the 'new' singularity provided here  works like charm.

Virulign does this for virus sequences, the publication is available here and github address here

Hi @kristoftheys welcome to the site, congratulations on Bioinfo paper. I've edited the post to include publication and github address.

@Michael many thanks for the completion. sorry for the inconvenience, next time better :)


text:
State of the art mutation simulation software
There are many features affecting mutation probabilities, e.g. CpG mutations are 10-fold more likely than other types of mutations.
Is there a model (preferably with software) which can take two aligned genomic regions, estimate parameters of the neutral mutation procces and then simulate mutations for another region? Obviously, there are models like HKY and GTR, but they are somewhat simplistic and do not take the context into account.
One model I found is this, but it seems that there is no software implementing the model. Is there something more recent?
I am also looking for software which can account for indels.

How good the simulation must be? Should it include that in transcription factors and translated genes are less mutated aside from CpG and indels transverstions...? (I am not aware of any software, but I think that with a single reference any tool will have problems estimating good parameters)

@Llopis the idea is to use it on all the regions with very weak selection (e.g. intergenic non-functional), only to account for mutations process, not selection, which of course affects genes and regulatory elements.

I'm not sure these tools allow parameter estimation from the data, but I've heard good things about Indelible https://academic.oup.com/mbe/article/26/8/1879/980884 and Evolver https://www.drive5.com/evolver/EvolverUserGuide.pdf

Does this help:
https://www.biorxiv.org/content/early/2017/11/22/223297 New synthetic-diploid benchmark for accurate variant calling evaluation

FoldX is the state of the art method for mutation modelling (cited by countless papers).  It does not do exactly what you want without extra work though.

The MMOS company developing Citizen science solutions for these kind of tasks.

Question I suspect this question is about the phylogenetics of methylation and the approach the investigator is proposing would be the last approach to use.
Summary
The approaches to assess the phylogenetics of methylation in order of preference are:

dN/dS between CpG sites and non-CpG sites,
An explicit molecular clock between CpG sites and non-CpG sites
A likelihood ration test (LRT) based on a null distribution of mutations generated by a Monte Carlo algorithm

Background/Rationale
The approaches you are using have been around for a long while and generate the random distribution of mutations for a given phylogeny. Its phylogenetics in reverse and the technique is known as Monte Carlo, you start with a randomised probability and send it through a parameterised model to predict the amino acid/nucleotide. Thus it is ML(max likelihood), Bayesian, "phylogenetic-HMM" in reverse. It is used within likelihood ratio test calculations. The model is determined by a standard ML, Bayesian phylogeny algorithm, i.e. its circular mathematics because there is no independent means of calculating the mutation behaviour, so the precise context you are using this for needs to be carefully considered.
Packages There are lots of packages that do this type of Monte Carlo simulation, SOWHat being a good one (LRT) and you generate a large number of replicate data sets (100 or 1000). The one of the base algorithms is "seq-gen", although PAML may have implemented this.
Considerations To use this style of approach you need to carefully consider your question. The "mutation simulations" when analysed via a phylogenetic program will produce the same parameters that you initially set and trees of basically the same length. If you are using this to generate a null distribution for a phylogenetics test these approaches are useful, you then compare the observed likelihood against the null distribution. If you are using it to work out whether the mutation rate at CpG sites is higher than at other sites, its one approach amongst a number of alternatives.
Drawbacks Calculating a de novo null distribution of the mutation rate is computationally very expensive and therefore tends to the last calculation performed. It will address a single hypothesis, usually regarding topology.
Crux of the problem To really do what (I think) you want to do you need an independent measure of mutation and mutational behaviour and its not trivial to achieve. You would have to consider it like a machine learning calculation, with a formal training set, I don't know whether such an approach has been implemented.
Strengths/Summary In summary, Monte Carlo phylogenetics simulation has a limited application because of circularity of the calculation, BUT, BUT, BUT what I have omitted is when it is appropriate it is a very powerful test indeed.

Points of interest I do like the population dynamic simulations within Beast, this certainly "a priori" and its potential IMO needs further exploration. However, I don't think you looking at how molecular epidemiology might impact the mutational behaviour methylation. I have not looked at Indelible, but looks interesting.
I was going to talk about "processor war simulations" but its probably off-topic.

Thanks for your response. My question is not so much about phylogenetics of methylation. What I need essentially is a null distribution for neutral mutations. MCMC or maximum likelihood helps a lot given that there is a model :) But I am not aware of a non-time revesible model accounting for sequence context and indels.

I don't know what "sequence context" is but non-reversibility of mutations (and by inference time) is Beast, all other ML or Bayes (MCMCMC) algorithms have a symmetrical mutation matrices. MCMC in Bayes is different from the Monte Carlo used for nucleotide simulation. To honest the lack of clarity in the question is of concern, so I can't help any further.

In classical phylogenetic models a probability of mutation at position N depends only on the nucleotide at that position and some measure of time. What I mean by sequence context here is the fact that probability of mutation will also depend on nucleotides and mutations at positions N-M..N+M. We have a wealth of empirical evidence that there is such a dependency. But not many models accounting for it (such as this one).
I will try to clarify the question and to give some background information there.

So do I, but I bet you can't name the mechanism.  Anyway Nick Goldman did this years ago, its not new


text:
What is the best way to account for GC-content shift while constructing nucleotide-based phylogenetic tree?
Let's say I want to construct a phylogenetic tree based on orthologous nucleotide sequences; I do not want to use protein sequences to have a better resolution. These species have different GC-content.
If we use a straightforward approach like maximum likelihood with JC69 or any other classical nucleotide model, conserved protein coding sequences of distant species with similar GC-content will artificially cluster together. This will happen because GC-content will mainly affect wobbling codon positions, and they will look similar on the nucleotide level.
What are possible ways to overcome this? I considered the following options so far:

Using protein sequence. This is possible of course, but we lose a lot of information on the short distance. Not applicable to non-coding sequences.
Recoding. In this approach C and T can be combined into a single pyrimidine state Y (G and A could be also combined in some implementations). This sounds interesting, but, first, we also lose some information here. Mathematical properties of the resulting process are not clear. As a result, this approach is not widely used.
Excluding third codon position from the analysis. Losing some short-distance information again. Also, not all synonymous substitution are specific to the third codon positions, so we still expect to have some bias. Not applicable to non-coding sequence.

It should be possible in theory to have a model which allows shifts in GC-content. This will be a non time-reversible Markov process. As far as I understand there are some computational difficulties estimating likelihood for such models.

I would just add that I think there's a key assumption in the setup here: "I do  not want to use protein sequences to have a better resolution". We can decompose 'better' here - it's likely to be more precise but also more biased, the latter for all the reasons you outline.

In case you might be interested, I tested some of the approaches you mention, plus a few other recoding schemes (http://dx.doi.org/10.6084/m9.figshare.732758) in the following papers: http://arxiv.org/abs/1307.1586 and http://dx.doi.org/10.1093/molbev/msu105

There are models that take into account compositional heterogeneity both under the maximum likelihood and Bayesian frameworks. Although the substitution process is not time-reversible, the computations are simplified by assuming that the instantaneous rate matrix can be decomposed into an "equilibrium frequency vector" (non-homogeneous) and a symmetric, constant exchange rate matrix.
I guess all your suggestions are also valid, and I remember recoding being used successfully to reduce the GC-content bias (examples in the references above and here).

The following 2004 paper describes a way to model compositional changes across the tree, in a Bayesian framework: https://doi.org/10.1080/10635150490445779
A python package implementing this ("p4"), and improvements added along the years, is available here: https://github.com/pgfoster/p4-phylogenetics
To get started, you may find useful examples here: http://p4.nhm.ac.uk/scripts.html
This has been used in a few large-scale phylogenetic analyses.

The answer is the logDet algorithm was constructed to overcome GC% clustering.
At the time it was devised only a distance method was available/implemented, so it wasn't very powerful. The posts here imply that a Bayesian or ML approach are available and these do tightly stick to the model.
Original publication here

Do you have a publication or webpage in mind? Can you link it?

Link provided above. Goes back a long way ... to 1996


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:
Testing for viral/bacterial infection FASTQ files
What's the best way to test if a pair of Illumina FASTQ files from blood DNA extraction contain any reads are from viral or bacterial infection from human samples? Thanks in advance.
EDIT:
Looking at the available options, I started testing Kraken, the metagenomics classification (Wood and Salzberg Gen Biol 2014) for samples where I have seen an indication of viral reads. I will update with results / further questions.

There are some people voting to close this question as too broad. Perhaps you could indicate what are your requirements and what options do you consider

I think you could try pathoscope for this.

If you're talking RNA, I have a blog post from a while back about Diamond/MEGAN:
http://achri.blogspot.com/2017/04/data-for-free-viral-metagenomics-in.html
In principle this will work for DNA if you have a non-trivial number of pathogen reads (since they tend to be gene dense in their genomes). Otherwise I use Kraken on a database for all the NCBI nt database, but this is slower and doesn't let you sanity check the false positive taxonomy assignments as well as I describe in the post. Of course, ideally, you run multiple methods if you have the resources to do so.

You may find this blog post helpful: Identifying viral sequences in TCGA data using Kraken and Centrifuge.

Our goal was to evaluate the performance of Centrifuge and Kraken for
  detecting viral transcripts across a cancer cohort. Centrifuge and
  Kraken produce very similar results and consistently detect HPV
  strains in cervical cancer.

There is also a review Overview of Virus Metagenomic Classification Methods and Their Biological Applications, but they focus specifically on viruses:

We evaluated 49 published computational classification workflows for
  virus metagenomics in a literature review. To this end, we described
  the methods of existing workflows by breaking them up into five
  general steps and assessed their ease-of-use and validation
  experiments. Performance scores of previous benchmarks were summarized
  and correlations between methods and performance were investigated. We
  indicate the potential suitability of the different workflows for (1)
  time-constrained diagnostics, (2) surveillance and outbreak source
  tracing, (3) detection of remote homologies (discovery), and (4)
  biodiversity studies.

To know if there are bacterial genomes in the sequenced sample I have used an approach based on kmers analysis. An useful tool to do that is kmerfinder which uses a kmers database built from bacterial genomes. Quoting from the kmerfinder web site:

The method is used to find the best match (species identification) to
  the reads in one or more fastq files or one fasta file in a (kmer)
  database produced using the KMA program.

I haven't use it for viruses genome analysis, but conceptually it should be work as well.
Hope this would be helpful.

What type of NGS data is it? If smallRNA or RNAseq data, you can simply assemble the reads into longer contigs and do the blast search using your reference viral or bacterial sequences.