text:
"perl: warning: Setting locale failed." in RepeatMasker
I'm trying to run Repeatmasker in Linux on the command line with:
./RepeatMasker -species human -alu -gff -dir /mnt/lustre/users/Analysis/RepeatMaskerOutput /mnt/lustre/users/Software/STAR/bin/Linux_x86_64/HRG_Newest_and_Annotation/GCA_000001405.28_GRCh38.p13_genomic.fna

but I am getting this error:
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_GB.utf8",
        LC_ALL = "",
        LC_CTYPE = "en_GB.utf8",
        LANG = "en_GB.utf8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker version open-4.0.9
TRF program not configured! This version of RepeatMasker requires a
local installatoin of TRF.  Please visit:
               http://tandem.bu.edu/trf/trf.html
to obtain the current version.  Once installed please re-run the
RepeatMasker configure script to setup RepeatMasker to use the new
installation. at ./RepeatMasker line 442.

I have downloaded rmBlast and the latest TRF file. I have put the TRF file in my Users/username/bin PATH folder and used chmod a+x (and chmod 755) on it. I then changed my LD_LIBRARY_PATH using export to my glibc directory.
After this I used perl ./configure and directed the program to the PATH folder where the TRF file was. It gave a series of the same error messages as above but was 'successful' as shown here:

[user@nodec23() RepeatMasker]$ perl ./configure
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker Configuration Program

**TRF PROGRAM**

  This is the full path to the TRF program.
 This is now used by RepeatMasker to mask simple repeats.

Enter path  [ /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/TRF.pm ]

Add a Search Engine:
   1. Crossmatch - Phrap             [ Configured ]
   2. RMBlast - Improved NCBI Blast  [ Configured, Default ]
   3. HMMER3.1 & DFAM                [ Un-configured ]
   4. ABBlast/WUBlast                [ Un-configured ]
   5. Done

Enter Selection: 5
 -- Setting perl interpreter...
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/repeatmasker.help

I'm at a loss of what to do next. I have googled and tried changing each of the locale settings with export to the following:

en_GB
en_US.utf8
en_GB:2000
en_GB.utf8

but none have changed the error.
I am running a Windows 10 machine and using PuTTY to connect to an HPC which uses a Linux environment.
Locale settings are as follows:
locale: Cannot set LC_CTYPE to default locale: No such file or directory
locale: Cannot set LC_MESSAGES to default locale: No such file or directory
locale: Cannot set LC_ALL to default locale: No such file or directory
LANG=en_GB
LC_CTYPE="en_GB"
LC_NUMERIC="en_GB"
LC_TIME="en_GB"
LC_COLLATE="en_GB"
LC_MONETARY="en_GB"
LC_MESSAGES="en_GB"
LC_PAPER="en_GB"
LC_NAME="en_GB"
LC_ADDRESS="en_GB"
LC_TELEPHONE="en_GB"
LC_MEASUREMENT="en_GB"
LC_IDENTIFICATION="en_GB"
LC_ALL=

Also please note that I do not have administrator access.
Typing 'perl -v' gives the following:
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

This is perl, v5.10.1 (*) built for x86_64-linux-thread-multi

Copyright 1987-2009, Larry Wall

Perl may be copied only under the terms of either the Artistic License or the
GNU General Public License, which may be found in the Perl 5 source kit.

Complete documentation for Perl, including FAQ lists, should be found on
this system using "man perl" or "perldoc perl".  If you have access to the
Internet, point your browser at http://www.perl.org/, the Perl Home Page.

Typing "LC_ALL=C perl ./configure" and going throught the configurations leads to the same outcome but without the perl errors:
Enter Selection: 5
 -- Setting perl interpreter...

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/k1894998/Software/RepeatMasker/RepeatMasker/repeatmasker.help

However, the error remains when trying to run Repeatmasker.

I think it is quite simple you need to install perlbrew here
You then run the software under perlbrew it should work.
At present you are running the system perl5 and that is causing all sorts issues, for example the age of the installation and the fact system operations and your software are using the same version of perl.
All the instructions are in the website link.

Actually I'd recommend Anaconda/conda now.


text:
Help me to calculate Heaps Alpha value from the roary pangenome pipeline result?
I need to know whether my pan-genome is open or closed. For that, I need to calculate the Heaps Alpha value from my roary pipeline results. Please help me to do the same.
Many thanks in advance.

Directly quoting from the link provide:
Roary By Andrew Page based on version 3.11.2 (22-Jan-2018)
Roary is a high speed stand alone pan genome pipeline, which takes annotated assemblies in GFF3 format (produced by Prokka (Seemann, 2014)) and calculates the pan genome. Using a standard desktop PC, it can analyse datasets with thousands of samples, something which is computationally infeasible with existing methods, without compromising the quality of the results.

I edited the question to give details of Roary. I personally find the context difficult to follow and further details would be helpful, such as the aim of the project as well as the genomes being used.

If this is a homework question (or otherwise canned question with a known answer), could you please tag as 'homework'? Also, please provide additional context around what you're trying to do (e.g. sample type, project goal).

source: Stack Exchange

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

sort: votes


text:
how can i download rna-seq data sample for starting rna-seq analysis by R?
I am looking for an induction into RNA-seq analysis using R. I have yet to decide a data set, or a  subject area. Any help in this regard would be welcome.

Hi @Aminafzalifar welcome to the forum! My advice is have a think about what to focus on and return with a better defined question. Can you code, do you understand dataframes? Sauret and Bioconductor are popular packages in R which deal with this, however in my opinion (not universally shared) you need a local supervisor to help you.

I'm sorry, but we don't do this sort of open-ended discussion here. This site isn't a forum, it's strictly a questions & answers site and we only cover specific questions that can have specific answers. Your question is just too broad for this site's format. You might want to take the [tour] or visit our [help] to understand how the site works. Especially this page that explains the kind of question to avoid.

thanks and sorry about my bad question.

You must have found the airway dataset by googling "RNASeq R tutorial"
https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html
Why won't this work for you?

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


text:
Global or local alignment on sequences for which it is assumed that they have common ancestry
As the title says I would like to know which alignment method(global or local) should I use for sequences for which is assumed that they have common ancestry?
The sequences are 21 aa long peptide chains (proteins).
About the background of the sequences, it's only stated that it's ASSUMED that they are from the same ancestry.
The thing that we assume that they are from the same ancestry makes me confused. I know that on similar sequences with the same length usually, we use global.

Depends. What question are you trying to answer? Please [edit] your question and give us some background. Also clarify what kind of sequences you have (nucleotides? protein? genes? cDNA? whole genomes? introns?). How large are your sequences? We can't help if you don't give us more details.

@terdon updated my question with more info. Its more of game with words.. the ancestry states that we should use global but the assume parts gives us idea of using local. Thats the confusing part.

Thanks for the edit, but you still haven't told us what question you are trying to answer. Why do you want to align them? What information do you want to get from the alignment? Are you trying to see if they are homologous? Are you looking for shared domains?

@terdon I am beginner in this field. I just want to know which is the right method for this. And yeah i would like to test if they really are homologous ( not just assume).

Thanks, I was asking because the method will depend entirely on the question you are asking :)

If you are trying to figure out the proteins' evolutionary history (e.g. if they are homologous) then you want a global alignment. In fact, there are very few cases where a local alignment is better, local is usually done only because it is far more efficient and less computationally expensive than global. But with tiny sequences of just as few residues such as the ones you have, you would want global almost every time.
The only case I can think of where you might prefer local alignment is when you have sequences that share little overall sequence similarity but you suspect they carry some sort of shared motif. For example, in proteins, that could be a protein domain. In that case, you would want to do a local alignment (possibly using something like blast) to identify those sub-regions of higher similarity.
Note that simple sequence similarity isn't really enough for you to decide whether they are indeed homologous. A more sophisticated approach such as psi-blast or HMMER will give better results. Nevertheless, a simple alignment will help you at least discard some obviously non-homologous sequences.

My initial reaction is that 21 amino acids is too little information to demonstrate these proteins are homologues (which is the question behind your question). A homologous site means it shares a common ancestry and is distinguishable from random chance or convergent evolution. A common ancestor can be demonstrated using phylogeny (trees).
If a regular alignment program cannot align these sequences, such as clustalo or muscle, the assumption they are homologues cannot be easily tested using phylogeny. That would then become complicated at which point you would need to ask whether it is worth pursing the question.
My recommendation is to perform a simple protein-protein blast and assess the e-scores to get some feel about you data. Your question is achievable, but you would solve this in context to additional sequence data. Whether that sequence data is a denser sample or else contiguous sequence data needs further exploration.


text:
Confusion related to effective population size
I am making a demographic model for African American population in which I say 80% of ancestry comes from African population and 20% comes from European population. I am using the time/population time from Gazave et al. (2014).
The authors doe not mention the effective population size (Ne) of African American population after admixture. Is there a publication someone can point me to that mentions the effective population size of African American (AA) or even European American (EA). I know there is not much difference between the two. It's just I need a figure to have my simulation up and running.
Note Ne is an estimate of the people size based genetic diversity and migration and is associated with the Fst statistics

I am a bit confused about what are you actually up to. I suppose that the admixture of African American and European Americal populations is a continuous process, what do you mean by effective population size after admixture?

That is a fantastic question and in my opinion a really important question in context. The Ne of Africa will be much bigger than the European population, however I'm less sure about the West African population which gave rise to the Afro-American population. I suspect it is still large because of the "out of Africa" theory.
My advice is to start by taking a complete look at the entire subject, from Tenesa et al (2007).  Ne is a really important parameter in any population study and extremely useful. I have feeling this could be complex, but I cannot hypotheses any further because I don't do human genetics (except for immunology).

Thank you Micheal for your answer.

authors: Антон Пушкин


text:
Is it possible to get such a plot
I wish to reproduce a very informative plot in google

Given a group of samples and a list of genes for which we have data how each genes has mutated and whether the alleles are homozygous of heterozygous. These are expressed as boolean values
Genes enriched in 2 pathways like cell cycle and p53 as shown in this figure
Map of functional alterations for a group of patients. Genes (rows) encoding components p53–DNA repair pathway;  are affected by selected functional events (percent of samples altered and types of alteration are represented by colored squares) across group of samples in column. Alterations of the pathway are observed stacked green bar plots at bottom
I have the boolean matrix of mutated genes per sample would be grateful for assistence.
An example of my copy number data
> dataset_gistic$annotations
     type                event
G2   "Homozygous Deletion"   "ACVR2A"
G13  "Homozygous Deletion"   "HIST1H3B"
G73  "Heterozygous Deletion" "CCNE1"
G74  "Heterozygous Deletion" "TSHZ3"
G78  "Low-level Amplification"    "CTNNB1"
G110 "Low-level Amplification"    "CCNE1"
G111 "Low-level Amplification"    "TSHZ3"
G119 "High-level Amplification"   "PIK3CA"
G121 "High-level Amplification"   "NIPBL"
G147 "High-level Amplification"   "CCNE1"
G148 "High-level Amplification"   "TSHZ3"
>

> dput(dataset_gistic$annotations)
structure(c("Homozygous Loss", "Homozygous Loss", "Homozygous Loss",
"Homozygous Loss", "Homozygous Loss", "Homozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "ACVR2A", "HIST1H3B", "CDKN2A", "PTEN", "NAV3",
"LIN7A", "ARID1A", "ACVR2A", "SCN3A", "CTNNB1", "PBRM1", "EPHA3",
"POLQ", "PIK3CA", "SLIT2", "NIPBL", "MSH3", "APC", "HIST1H3B",
"ARID1B", "CDKN2A", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "NAV3", "LIN7A", "SIN3A", "TP53", "ERBB2", "GPATCH8",
"RNF43", "GATA6", "SMAD4", "CCDC102B", "STK11", "SMARCA4", "CCNE1",
"TSHZ3", "CTNNB1", "PIK3CA", "SLIT2", "NIPBL", "HIST1H3B", "ARID1B",
"EGFR", "KCNQ3", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "SIN3A", "ERBB2", "GPATCH8", "RNF43", "GATA6", "CCDC102B",
"STK11", "SMARCA4", "CCNE1", "TSHZ3", "PIK3CA", "NIPBL", "EGFR",
"KCNQ3", "CCND1", "KRAS", "ERBB2", "GPATCH8", "RNF43", "GATA6",
"CCNE1", "TSHZ3"), .Dim = c(77L, 2L), .Dimnames = list(c("G2",
"G13", "G17", "G19", "G24", "G25", "G38", "G39", "G40", "G41",
"G42", "G43", "G44", "G45", "G46", "G47", "G48", "G49", "G50",
"G51", "G54", "G55", "G56", "G57", "G58", "G59", "G60", "G61",
"G62", "G63", "G64", "G65", "G66", "G67", "G68", "G69", "G70",
"G71", "G72", "G73", "G74", "G78", "G82", "G83", "G84", "G87",
"G88", "G89", "G90", "G92", "G93", "G94", "G95", "G96", "G97",
"G100", "G102", "G103", "G104", "G105", "G107", "G108", "G109",
"G110", "G111", "G119", "G121", "G126", "G127", "G131", "G132",
"G139", "G140", "G141", "G142", "G147", "G148"), c("type", "event"
)))
>

And this is sample of my mutational data the number of mutations for each gene in each sample
    om$numericMatrix
       LP6005500-DNA_D03 LP6008336-DNA_G01 LP6008460-DNA_D01 LP6005334-DNA_H01
TP53                   2                 2                 2                 2
CDKN2A                 3                 1                 2                 0
       LP6007600 LP6008334-DNA_A03 LP6008334-DNA_A04 LP6008334-DNA_B02
TP53           2                 2                 2                 2
CDKN2A         0                 0                 0                 0
       LP6008334-DNA_C02 LP6008334-DNA_D02 LP6008336-DNA_H01 LP6008337-DNA_A07
TP53                   2                 1                 2                 2
CDKN2A                 0                 0                 0                 0
       LP6008337-DNA_H06 LP6008460-DNA_A04
TP53                   2                 2
CDKN2A                 0                 0
>

> dput(om$numericMatrix)
structure(c(2, 3, 2, 1, 2, 2, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2,
0, 1, 0, 2, 0, 2, 0, 2, 0, 2, 0), .Dim = c(2L, 14L), .Dimnames = list(
    c("TP53", "CDKN2A"), c("LP6005500-DNA_D03", "LP6008336-DNA_G01",
    "LP6008460-DNA_D01", "LP6005334-DNA_H01", "LP6007600", "LP6008334-DNA_A03",
    "LP6008334-DNA_A04", "LP6008334-DNA_B02", "LP6008334-DNA_C02",
    "LP6008334-DNA_D02", "LP6008336-DNA_H01", "LP6008337-DNA_A07",
    "LP6008337-DNA_H06", "LP6008460-DNA_A04")))
>

Why unvoting? Do you think only by a simple googling one can get such a plot?

Can you update the question by giving an example of your data looks like instead of describing them by words?

@Exhausted, please provide the community with a reproducible example: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example As was suggested in another question of yours, dput() will help us "reproduce" your input data. If you invest a couple of minutes going through the link and subsequently increase the "quality" and "completeness" of your questions, you will get answers much faster you do now.

While it was easy to get your data after your edit, the data was not quite useful, the second data frame/matrix has 2 columns only!
I had to create a fake data set, that looks like this:
my_df <- data.frame(
  gene_class = rep(letters[1:4],3),
  some_value = rep(1:3*10,4)
)
rownames(my_df) <- paste("Gene-",rep(LETTERS[1:4],each = 3),1:4,sep="")

> my_df
        gene_class some_value
Gene-A1          a         10
Gene-A2          b         20
Gene-A3          c         30
Gene-B4          d         10
Gene-B1          a         20
Gene-B2          b         30
Gene-C3          c         10
Gene-C4          d         20
Gene-C1          a         30
Gene-D2          b         10
Gene-D3          c         20
Gene-D4          d         30

Using the ComplexHeatmap package, I got something similar to the figure you have posted.
# data frame for annotation
my_annot <- my_df[,"gene_class",drop=FALSE]

# the Heatmap() function requires a matrix
my_matrix <- as.matrix(my_df[,"some_value",drop=FALSE])

dev.off() # clear the plotting area in case not empty
ha = rowAnnotation(gene_class = my_annot$gene_class)

ht <-  Heatmap(my_matrix,
               left_annotation = ha, # gene_class
               row_names_side = "left", # position of row (gene) names
               heatmap_width = unit(0.3, "cm")*nrow(my_matrix), # cell width
               cluster_rows = FALSE, # no clustering
               cluster_columns = FALSE, # no clustering
               cell_fun = function(j, i, x, y, w, h, col) { # add text ("some value column") to each grid
                 grid.text(my_matrix[i, j], x, y)
               })
draw(ht, heatmap_legend_side = "bottom")

The ComplexHeatmap package offers enormous flexibility for customization. I strongly recommend devoting some time to it. The way I use it, I try to put all of my data/information in a single data frame/table, which I subsequently subset to create data frames and matrices for annotations and the actual heatmap(s) respectively.

Thank you @haci . You are amazing and helpful as always, definitely your solution works but I don't know how to format my own data to work with your code. What you have produced is what I need as long as I am able to adapt my data our your dummy input


text:
Alternative splicing positions on protein
I have mass-spectrometry data (Human cellsanalyzed and annotated by MaxQuant).
I want to find which peptides are covering positions of former splice junctions (on the mRNA sequence), based on protein id and peptide positions.
Where could I find an appropriate data? I want to know, based on peptide-sequence, annotaion and position on protein, whether this peptide is a single exon, or if it contained areas of former splice junctions.
Itai

Hello Itai, what do you mean by "appropriate data'? Do you mean where you can find annotations of splice junctions? What species? Please, [edit] your question and be clear and verbose. The second thing that could be improved is that your question shows very little research effort. What have you tried? Have you found any clue how the problem could be resolved? Yes? Write it here.

I edited the question

Well done, It's way clearer now. Maybe just add the species you are looking at and the results of your own research on the topic...

If you are considering the human genome as "data" (not mRNA), I would run your single peptide via tblastn which will take your protein(s) translate them to nucleotides and run a blast against Genbank. The splice sites will be immediately apparent because they force indels in your query. If there are no indels there are no introns (splice sites). If you want to process several proteins, one solution is to use eblast, a command line version of tblastn  (Anaconda package). You can confirm this against a full genome annotation (note, I don't work on humans).
Alot of work has been done on splice sites, a friend worked on this, but I think you just want to focus on a restricted number of target peptides in your mass spec. output.
mRNA should be avoided IMO if you are looking for splice sites, i.e. intron boundaries.

eblast sounds like a great idea, thanks!


text:
What do negative dN or dS mean in cal_dn_ds in Biopython?
I'm noobie in bioinformatics. I use cal_dn_ds from Bio.codonalign.codonseq with method 'NG86'. Sometimes I get dS=-1. What does it mean? How can number of substitutions be negative? Does it mean that something is wrong with my data?
Upd.
1) Chris_Rands pointed to sources, I read them already, but the question remain. What is that ps variable, and why if it's greater than 3/4, dS = -1?
2) Some more details.
cal_dn_ds requires sequences from CodonAlignment object, which is created by function Bio.codonalign.build that requires (quote from docs) protein alignment and corresponding nucleotides. I have nucleotide sequences from NCBI at the beginning. So I trim non-coding regions from nucleotide sequences, translate all trimmed sequences, write them to a file, perform multiple alignment of the protein sequences, and then pass the list of aligned protein sequences and the list of nucleotide sequences to the build function.
I analize protein DKK1, and example sequences that get negatives are Empidonax traillii XM_027898752.1 and Physeter catodon XM_007110600.1. But I just want to understand general idea behind dS = -1, what can it mean it different cases.
3) Michael G. gave an interesting suggestion that it means zero, but it's strange because the function returns dN and dS separately, so Biopython might as well return zero, so the user would handle those cases himself. And for some sequences Biopython actually returns normal positive dS and dN = 0 at the same time. If it's an error, what are the possible causes? And why in those error cases Biopython can calculate dN, but dS is -1?

Please [edit] your question and give us some more details. What kind of input do you have? What are you passing as input to cal_dn_ds? If you are using some sort of alignment, what species are involved? The more details you give us, the better your chances of getting a useful answer.

https://github.com/biopython/biopython/blob/master/Bio/codonalign/codonseq.py#L390

terdon, Chirs_Rands, thank you for replies, I updated my question with more details.

dS = -1
.. there are two suggestions

I suggest its an error, possibly based on nucleotide calls within the S, check for anomylous nucleotides, particularly over a non-synonymous codon
It might mean their are no synonymous mutations. It is not the best idea to assign this a value of 0 in a ratio because it will cause the ratio to fall over when in fact their could be strong selection pressure, ie non-synon mutations. Having an incrementally small value is problematic because dN if positive could anything up to infinity, so -1 might flag this problem.

Okay I now know why... You are correct to be concerned with the result because its serious.
My sincere advice to you is do not use any dN/dS ratio where -1 is occurring under any circumstances. On this occasion I'm not going to explain why for two reasons:

First and foremost, as a phylogeneticist you should know why, it relates to the fundamentals of the discipline.
(Hint) The answer is in my past posts on this forum.

To conclude, I'm pretty nuanced in statistical dogma, you can check my record, IMO most dogma is tribalism because every calculation has weaknesses. HERE HOWEVER, there is a real problem and again I urge you to ditch any data set where -1 is occurring and everything associated with it. That is a fairly dire warning should give you a clue about whats wrong. If you are not prepared to do this, then you must report the result and the situations where it is arising (don't hide it).

Thank you for your reply, I provided some more details for you at Upd 3. I checked my nucleotide sequences, and they do not have anything except ATGC, if that's what you meant by 1.


text:
COVID-19 GWAS: what is a cross-replicating association?
In a genomic study of patients infected with SARS-CoV-2, the authors  detected cross-replicating associations with rs11385942 at locus 3p21.31 , the association signal spanned the genes SLC6A20, LZTFL1, CCR9, FYCO1, CXCR6 and XCR1
I would like to know what's a cross-replication association.
Source pdf is here

Replication in GWAS is an analytical term, approximating to 'robustness'. It is not biological phenomenon of a genome replicating. What the authors are saying is these genes associate with an allelic variant(s) which correlate with the disease phenotypic of COVID-19.
Three of those genes are very interesting indeed notably CCR9, CXCR6 and XCR1, because these are chemokine receptors. The later for example is associated with regulating between the innate and adaptive (CD8+, i.e. cytotoxic T-cell) immune responses (described in the image below). The innate/adaptive interface is an extremely trendy area in immunology.
The authors are implying the immunological signalling associated with one or more of these chemokines is fundamental to combating COVID-19 infection. 'Replication' is simply the GWA geek for 'its real'.


text:
Why is there a chloride ion in this 3D model?
Here we have a 3D model of a sequence that match with this sequence
VVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIK
in Blast
the 3D structure is here
Why is there a chloride ion in this 3D model?

It is a crystallisation artefact. Namely, the protein are placed in a condition where they fall out of solution without aggregating a gloop (as happens in most well in xstal trial) and to pack as a crystal they need to be placed nice and orderly. It is coordinated by water so is not relevant. A lot of structures have these —DMSO and ions are the most commonly.
Two footnotes.

1ZV7 is a crystal structure, i.e. solved by X-ray diffraction. A model is generally used to described an in silico model, such as those from Swissmodel, phyre, I-Tasser etc. These are generated ab initio (from scratch) or threaded (mapped on top of) against crystal/EM/NMR structures. For Covid I-Tasser has a nice picture of its models.
The best source of info about structures is the PDB. RCSB PDB is more commonly used, while PDBe is more tricker to navigate (or used to be). Both have the same data, but are better for the extras than NCBI.


text:
NCBI Blast databases - cannot find some databases
I'm working with a system that does local BLAST searches using some preconfigured BLAST databases.  From what I see, these were updated within the last year, but I cannot find any current databases available that have those same names.  One specific set is a non-redundant set of human protein or nucleotide sequences.  The databases are named human_nr and human_nt (with various extensions as per the BLAST db format).  There are others names up_sp_human (Uniprot-Swissprot human?) and up_tbl_human (Uniprot-Tremble human?), but a cannot find any similar databases on the NCBI site or other sites.
Where else should I look for such preconfigured databases?

I don't think there is a human_nr or human_nt database for NCBI.
An equivalent database is accessible via the standard nr/nt database and simply restrict the organism parameter to human or to taxa id 9606 . This essentially will do the same thing.
Or maybe this sort of thing?,
blastdbcmd -db nr -entry all -outfmt "%g %T" | awk ' { if ($2 == 9606) { print $1 } } ' > human_flat-list.db

This is super helpful!  I'm guessing these were custom prepared similar to what you describe.  Thank you!

Thanks for your feedback @KirkD-CO


text:
lower mapping rates in salmon v0.13 compared to previous versions
Hi there :) Thanks for the tool!
I recently updated to the new salmon (from 0.8... its been a couple years) and I noticed that my mapping percentages change dramatically between the two versions.
For example, using the default settings in v0.8, I see a mapping rate of 96.25% however upon upgrading to v0.13, these mapping rates go down to 82-84% depending on which of the parameters I play with (I've permutated  the validateMappings, incompatPrior, consensusSlack, minScoreFraction, ma, etc). Below I have a snippet of one of these tests.  I'm currently assembling the unmapped reads to see if I can determine why these reads were not mapped.
What do you think is the most likely cause for this dramatic change in performance?
Perhaps it could be the transcript index? (I used a kmer of 31 and my reads are 100 bp).
UPDATE:
I have run the indexing and quasi-mapping  using salmon v0.8, v0.10 (when validateMappings was released), v0.12, v0.13 and had mapping rates of 96, 90, 90 and 83 respectively.  So it looks like a gradual decay of mapping rates.  Or alternatively, a gradual reduction in spurious mapping?
Seems like what ever improvements from 0.12 to 0.13 changed the mapping rates significantly.  I will go ahead with both mappings (v0.12 and v0.13) to see if there are differences in the DESEQ2 outputs.
Thanks a lot!
(while I have presented values here for one set of fastq files, the same pattern is seen throughout all 72 of my samples from different experiments)
for F in $FILES ; do
  R1=${F}_L001_R1_001.fastq.gz
  R2=${F}_L001_R2_001.fastq.gz
  salmon quant \
      -p 35 \
      -i $ref_dir/$transcript_index \
      -l A \
      -1 $R1 \
      -2 $R2 \
      -o $(basename $F).VM.rFB4.IP9.quant \
      --incompatPrior 9.9999999999999995e-21 \  # default of v0.8
      --writeUnmappedNames \
      --validateMappings \
      --rangeFactorizationBins 4

done

Update salmon is now released with a --allowDovetail option in  0.13.1 and inspected via a -z dump.

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed.

[Most likely answer based on comments attached to the question]
Hi Courtney, We've looked into the data and identified the source of the different mapping rates. Specifically, the cause is discordant reads; in this case reads that are dovetailing with respect to their positions on the reference. If you run Bowtie2 with "RNA-seq" flags (i.e. disallowing discordant mappings), then you get essentially the same mapping rate as salmon here. Is allowing dovetailing reads something that should be very important in your analysis? If so, we can consider what would be necessary to add back support, but generally, they provide evidence of tenuous quality.


text:
PAML branch-site model: Identical lnL value between null and alternative models?
I have been running PAML, via codeml, to try to identify positive selection in certain branches of my tree for one gene (but will eventually be trying to automate the analysis over multiple genes). This was done by setting model = 2 and NSsites = 2 in the codeml.ctl file. I attempted to run both the null and alternative models by setting fix_omega = 1 in the null.ctl and fix_omega = 0 in the alternative.ctl.
However, after running both models, they have returned exactly identical log-likelihood (lnL) values. I wasn't expecting to find anything significant, but the fact that they were identical struck me as odd. Out of curiosity, I altered the branch of interest in the newick tree file and reattempted the analysis, but again, the lnL values remain the same between the null and alternative models (although slightly different to the value returned with the different branch selected).
I have since repeated this many times, choosing a different branch(es) each time to see if I can get a difference between the null and alternative, but no luck. Am I missing something obvious? Should I really be getting the same exact result between null and alternative models, regardless of where in the tree I place #1? To reiterate, I wasn't expecting to find significant results, but was at least expecting (non-significantly) different numbers for lnL...

The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.
Possible explanation fix_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.
Goto debugging As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,
time codeml ./file

The thing to remember is fix_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.
Nitty gritty debugging
1 Please check you have
seqtype = 3

If not, then that is your answer, omega was never being set.
2 For debugging only make sure you have,
model = 1

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below
3, Then try a different NSsites model,
NSsites = 8

This should produce real differences. The NSsites options are,

0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w9:betaγ
10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1; * 13:3normal>0

w here is for omega
4 The proof of your system is that when you
fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

This should produce likelihood variation

Finally
model = 2

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it.  I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.

Thank you Michael, you've given me something to try. If I'm still struggling after debugging, I'll update this question.

Any updates @PollardMD ?

Apologies for the late reply. I managed to get it working eventually. It looks like some genes/branches that I am looking at truly do produce identical values for both null and alternative models. When looking at other genes of interest, I have been able to get non-identical results from the null and alternative models. I guess I just wasn't expecting the lnL values to be perfectly identical, and that threw me off.


text:
I am getting a parsing error when uploading .chip file to GSEA software
I get the following error when I try to load a .chip file to GSEA software. What is the problem?
---- Full Error Message ----
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "YBL028C"------
java.lang.NumberFormatException: For input string: "YBL028C"
    at sun.misc.FloatingDecimal.readJavaFormatString(Unknown Source)
    at sun.misc.FloatingDecimal.parseFloat(Unknown Source)
    at java.lang.Float.parseFloat(Unknown Source)
    at edu.mit.broad.genome.parsers.TxtDatasetParser._parseNoDesc(TxtDatasetParser.java:169)
    at edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:131)
    at edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:87)
    at edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:202)
    at edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:749)
    at edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:726)
    at edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52)
    at javax.swing.SwingWorker$1.call(Unknown Source)
    at java.util.concurrent.FutureTask.run(Unknown Source)
    at javax.swing.SwingWorker.run(Unknown Source)
    at java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
    at java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.lang.Thread.run(Unknown Source)

Edit:

I am also using GSEA 3.0. What does not make sense is the other files have the right extensions and are not getting this error, but only this file is failing to upload.

https://www.gsea-msigdb.org/gsea/index.jsp
"Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically
significant, concordant differences between two biological states
(e.g. phenotypes). "

Hi there, would you mind sharing a bit more details? What have you tried to resolve your problem? DO you think it's a problem of installation or the file? Note also that you can [edit] your question.

Also perhaps it could be useful to know which version of GSEA are using and how are the first lines from the file (or those around YBL028C )

Please show us an example of the file you are trying to load. It is complaining about YBL028C. That looks like a gene/protein name, but the tool is trying to read it as a number. So where do you have YBL028C in your data? What is a .chip file?

did you already solve the problem? I also came across the same error when uploading a similar .chip file.


text:
Restriction enzymes with multiple recognition sequences
In a bioinformatics programming course I am currently teaching (using Python, and approached from a mostly CS rather than a biology background) I assigned a project which involves finding sites where restriction enzymes bind to a DNA sequence. The project is heavily influenced by Chapter 9 of the now classic text "Beginning Perl for Bioinformatics" by James Tisdall. I am using the REBASE database, in particular REBASE Format #9. In this format, each restriction enzyme has a recognition sequence which looks like e.g. C^YCGRG where ^ gives information about the exact cut-point and upper-case letters such as Y and R other than ACGT are ambiguity codes. The point of the project is to turn these into regular expressions where cleavage information is discarded and the ambiguity codes are expanded into character classes, and then use these regular expressions to search for matches.
All fairly straightforward, but one very observant student noticed that some enzymes had multiple sequences, for example M.BssHIII has a total of 5 sequences: ACGCGT,CCGCGG,RGCGCY,RCCGGY,GCGCGC. It is easy to see that these correspond to a total of 11 sequences (after the Rs and Ys are expanded), far fewer than the 144 sequences which would match VSSSSB (which is the least ambiguous single sequence capable of matching the above).
My question is, what is going on here? Am I safe in simply replacing , by | (for "or") in the regular expression, or is there something more subtle going on? Only 8 of the nearly 5000 enzymes have such multiple sequences, and neither the Tisdall book nor even the description of the format (both of which seem to assume single sequences) shed any light on the situation.

It's programming exercise. Just use "|" for regex. It would be also biologically correct.

I would honestly just remove the non-homologous cutting site enzymes from the database. Firstly, these are biology students and require a simple picture to grasp the concept of coding via regex pattern matching. Secondly, in a wet-lab nobody uses these enzymes, I was not even aware the "8" existed until your post. In a teaching lab the best motto is "keep it simple".

I was tempted to remove these 8 from the database, and might end up doing so, but I was genuinely impressed with the student who noticed it and told her that I would look into it. Most of the students are actually CS students taking it as an elective, but I have a couple of biology majors or students that have some major-minor combo between CS and biology or chemistry.

I guess the cat is out of the bag. The difference is this is teaching rather than research. In my experience I would create a subproject to avoid distraction, and offer to all students is interested. Thus if you introduce the | and use it to filter the "anomalous 8 restriction enzymes" into a dictionary, panda, or np, to be analysed later if there is interest. Its interesting a core-bioinformatics text has a "bug".

At first I thought that these enzymes postdated the Tisdall book (which was published in 2001) but a check of the references for these shows that some of them were studied as early as the 1980s. I'll have to check out the book again, it is possible that there was something I missed. It is also possible that in an earlier edition of REBASE these were either missing (since when Tisdall wrote there was only around 1000 enzymes included) or had ? for their recognition sequence. I'll see if it is possible to download earlier versions of it. It could be interesting to see how it changed over time.

Keep in mind 1980s the way cleavage sites were determined different than today. Currently, it's easier to find degeneracy in enzyme cleavage site by NGS. You could also email Dr. Rich Roberts (Nobel Laureate)  and Dana Mecelis. I am sure they will get back to you.


text:
Calculate the percentage of each unique phylogenetic tree in a BEAST output
I have a nexus formatted BEAST output containing 20,000 phylogenetic trees of seven taxa. Is there any way to get the percentage of each unique phylogenetic tree contained in this output?
I already made an unsuccessful attempt with R.

If the phylogenetic trees can be represented as strings then it is basically just string matching. Anyway, I think your question would be more suited for [bioinformatics.se] as it is mostly about the programming aspect.

Could you explain what is a "unique phylogenetic tree"? Perhaps if you upload an image it would be clearer what do you mean.

@llrs When I mention a 'unique phylogenetic tree', I am referring to the topology of the tree. I would like to extract all the unique topologies from my BEAST output, and calculate how many times each of these topologies appears in my BEAST output. This would allow me to establish a percentage for each topology. The idea is to know which topologies are the most represented in addition to the consensus tree.

I finally managed to do it in R. Here is my code:
install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)

This is fantastics, thanks for posting back

You require a consensus phylogeny.
This is available in RAxML, PAUP, the old Phylip package and possibly the stats companion to Beast, Tracer. I'd need to dig up the code for RAxML and PAUP. The easiest is Consense here because you just download it and input the sequence. NOTE their is possibly a format difference (RAxML and Consense use Phylip format), I suspect you've got nexus, or Newick format.
Beast is likely to have its own consensus formating OR the Beast companion  "Tracer" may have this (I can't remember). Tracer is here and would be cool if it did consensus phylogeny because there will be no conversion problems (it doesn't use Phylip format). The alternative would be to use PAUP (no conversion problems), which is currently in a state of flux from various stages of freeware and updates. PAUP is contained within the Genious package, which is available as a 1 month free-ware.
Have a look and let me know: because if you want to use RAxML or Consense you'll likely need to remind me of the Beast output, if it is what I think it is a format conversion app is needed.  The calculation is easy, its easy to script a conversion program, but I think tree formating will be an issue.
You could use Ape, within R, but as you have gathered there are alot of other solutions. There are more recent, impressive, solutions as stand-alones which I've not used yet.

Many thanks for your answer! However, I already have a consensus tree and I am now trying to get the other most represented topologies from my BEAST output.

Ok. As far as I'm aware that program has not been written, at least a couple of years ago it wasn't, not sure about this now. Using R "ape" is your best bet. On this question I definitely can't help.


text:
Help interpreting my phylogenetic tree construction of bacterial species in same genus
I constructed a phylogenetic tree using mole blast, and it created a neighbor joining tree, so I used the alignment of those sequences with mine (1-43) in mega with TN93+G+I model but I don't understand how to interpret why my species in question (29) is so much different (or maybe novel?).
Attached is my original tree (Picture 1) and my bootstrap analysis (Picture 2):

Are these the same tree with different roots? The first tree is probably rooted by middle point. How did you root the second tree?

Yes that is correct. The interpretation is easy and relates to the population structure of two definite bacterial species within the same genus. I will provide a detailed interpretation shortly.

The top phylogeny is likely mid-point rooted, but I suspect the rooting is fine. The lower diagram has been rooted on the query sequence. It is taxonomic analysis using 16S which in bacteria is considered to represent species delineation at a given percentage of divergence. I thought it was 5%, which is impossible based on the phylogram in the top picture.
Bootstrap values lower than 80% mean that the clade is not supported (monophyletic group is a better term). The delineation between 10001 and 10013 against the massive ingroup of 13 16S sequences is robustly supported by bootstraps of 95% and 99%. The divergence between them is around 2%. There is no taxonomic distinction between the 13 16S sequences based on the phylogram.
One question is whether this is full length 16S? If it is there appears to be a species designation confusion in the lower clade supporting 10001 and 10013 because it is identical to to Dermococcus species and there is no way to assess which one is valid, nor whether they are really members of the same species. The very small divergence suggesting in D. abyssis might also call this into question. 10013 is in fact a robust outgroup to 10001 and the three Dermococcus spp. but the genetic distance involved is small.
However, this is not the major result. You have recovered what appears to be a large population (if that is how you sampled it) with close genetic homology to Dermococcus strain DSM20448. Given the complexity of species lineation within this genus it is impossible to say whether this represents the same species as DSM20448 but you could call it a "species complex".
Overall the population structure (if that is how it was sampled) is heavily skewed towards a bacterial species which is hitherto received minimal attention. This population skew (if that is how it was sampled), could be reversed in different ecotopes/niches/ecoregions/patients?
In summary, I wouldn't get caught up in systematics because in this genus its a bit messy. I would focus on the population structure and these results are very clear.

Thank you so much! It is only a partial sequence so the bootstrap being low isn't a surprise.


text:
How are the scores of GeneSplicer, MaxEntScan and SpliceRegion interpreted from VEP annotations?
I am using VEP from Ensembl to annotate my VCF files with the extra plugins of GeneSplicer, MaxEntScan and SpliceRegion. However, I don't fully understand the output of these scores.
I know that they are used for splicing predictions on trained sets.

What are the range of these values ?
What does a high or low value mean ?
What is the threshold ?

Can someone please explain them please ?


text:
How can `limma` be used to determine genes predictive of a given phenotype?
Using limma within Bioconductor here, which is a linear model package for microarray data and lmFit on gene expression data, I can perform a series of regressions gene ~ covariates by specifying design = model.matrix(~ covariates) and calling lmFit(eset, design).
How can I instead run the regressions phenotype ~ gene + covariates?

I think you can ask this question better at bioc support site. There the developers of limma can answer you.


text:
Should I compensate for two clusters of samples in microarray by hierarchical clustering and MDS?
I seem to have two groups of samples based on hierarchical clustering as well as MDS plot analysis using the sample label from the tree cut of the hierarchial clustering. I am analyzing microarray data. Should I adjust for the two different groups of samples in the model? For context, my covariate of interest is a continuous measure (air pollutant measure). This is a human cohort and not a case-control model and hence I was wondering if I should adjust for the different groups of samples.
[
MDS plot is here for the same samples:

Very interesting, if I understand correctly you are asking whether a straight hierarchical clustering should account for a continuous variable because it is resulting in discrete groups (clusters) which seems contradictory for a continuous variable?
If that is the question the answer is strictly no. Do not use weighting for this particular analysis (other regression style analysis yes - but not this one). Thus 'weights' are for regression analysis not for multivariate statistics. This clustering analysis is about assessing the correlation and not trying to impose it - thats regression.
The 'weighting' that can be used is different for a tree based analysis - it accounts for rate variation within the expression data - you are assuming rate homogeneity. That doesn't always work - but it's complicated and your data looks good as is.
Moreover, I would suggest your have 4 clusters. What you want to do is map the continuous trait across all branches 'tips' of your clustering. This will break up your continuous variable (hopefully) and looking at the MDS very much into discrete groups. That is the result, in fact it may not be 2 clusters but 4 clusters - map all the 'tips' against the continuous variable. There're statistically randomisation tree methods to assess the robustness of the observation - you can also use bootstrapping.
The question is whether these discrete groups are in fact the best description of the data. The MDS suggests that the continuous variable "air pollution" is in continuous within gene expression - using a combination of clustering and MDS you highlight this facet of the data. There alternative analyses namely regression or GLM to assess the relationship between the continuous variable and your results - the MDS is suggestive of that. However, it's not a 'and', 'or', 'right', 'wrong' relationship a continuous variable can be broken up into discrete groups: if there is a 1:1 relationship between the gene expression and the continuous variable "air pollution" (MDS analysis) - that is a better description, if that relationship doesn't hold breaking into discrete clusters is the best description thus far.
Essentially this is statistical modelling and the question is what is the best model to fit all the data. I'm personally not a fan of hierarchical clustering as a tree description of the data - but in your data it looks good.
At the heart of the question you are asking whether there are boundary effects in the data as described by the hierarchical clustering? Thus is there a certain level of pollution (or whatever) causing one gene expression boundary and then another and so on ... Perfectly possible - but central to that hypothesis is whether clustering is the best description of the data rather than a regression analysis, i.e. is it continuous or is it discrete.
Discrete vs. continuous  If it is discrete then k-means clustering would be a really, really good follow-up analysis to confirm the discrete "boundary" effects within the data. GLM is an analysis that might trap the data if it's continuous, there a lot of other methods in regression analysis. Basically if you map each "pollution" exposure against the MDS and if it's a perfect fit - trapping the result in a regression model outside eigenvalues would be the traditional thing to do. However, these days "unsupervised learning" (of which MDS is one method) is trendy so it may not be needed.
k-means clustering I would certainly consider k-means clustering in any case - this is because hierarchical clustering makes a rate assumption which might smash up the discrete groups. K-means can capture that without getting into complicated tree analysis. Basically, k-means you'll set 2 groups - assess the fit of the continuous variable, try 3 and then 4, how does the air-pollution variable fit? The data will also work well with supervised machine learning. However, it is genuinely continuous then regression analysis is best approach. With k-means you can also do parameterisation, i.e. whats the optimal number of groups represented by the data: if it comes back with '4' - which neatly segregate "air pollution" thats really good.
It's an exciting result. The way you write the result in the context of the conclusion is really important - i.e. the continuous verse discrete best fit analysis - in case you encounter a statistical modeller in peer-review.
My answer is a bit waffly, but is a representation of statistical modelling in this particular data, because two approaches may both be appropriate, but one is a better description than the other.

Note, if I've misunderstood the question and some sort of standardisation is needed between experiments - thats a different answer. Thus if you are saying that the groups represent different experiments, e.g. performed on different days or different arrays and there's risk of experimental artefact: thats about standardisation and the answer would be absolutely. Standardisation is essential to rule out experimental bias - particularly in microarray data.

Thank you for your very detailed answer. Highly appreciate it. I will slowly go through it and digest it. :)

Thanks let me know @siddharthadas


text:
Phylogenetics software that can represent collapsed nodes as triangles
I am trying to represent a large (100-300+ sequences) phylogenetic tree in a way that gets my point across. I figured that reducing the number of branches (by collapsing clades which include highly similar sequences) would be a good idea. The software that I'm currently using (iTOL) technically has an option for that ('collapse node'), but I'm not exactly a fan of the particular graphical representation, so I'd like to have a look at other options. Thank you for your time!

Try Figtree

For Mac Figtree here is highly recommended, it is also available for Windows and Linux here. The technique is called "collapsing clades" and is the only package I know that does it. I would recommend using it, because I would assume it is averaging the branch length of a given clade and the 'base' of the triangle is proportional to the number of taxa. Thus it recognised as being very reliable, albeit the precise algorhithm it uses ... I don't know :-O
With FigTree you get two different "triangle options" in the menu, which give different shaped triangles for collapsing the number of lineages within the select clade you chose. The big one, is the algorhithm appropriately described above. It has many options for producing publication quality diagrams, albeit I then import the output into a graphics program for further editing (there are shortcomings IMO in colouring and fonts. Then again for free is better than good.
If you are outputing a tree within a script, ETE3 (Python) will do a "collapsing of clade" option, but it isn't pretty triangle stuff. ETE3 will produce a publication quality tree diagram from this and reduce to 'low resolution', which is very useful in large analyses of loads of trees.
In summary, I don't think FigTree has any competition in the context you are seeking.
Moreover, you can configure your bin, or alias file so you can automatically open a treefile just typing,
figtree treefile.tre

Also you can reduce you font size to next point 1 or 2 (okay you can do that in anything), if you are attempting optical analysis on loads of taxa. Its also really stable. I just find it really useful, like a "jazz standard" its a "phylogeny standard".


text:
Should I expect / use a prokaryote-specific substitution model when building a phylogenetic tree?
I have a couple of prokaryotic genes that I'm trying to build phylogenetic trees from.
I inputted them into IQTree, which, as part of the workflow, runs through multiple substitution models until it finds the one that seems to fit the data best.
My issue: quite a few of the genes ended up with eukaryotic_specific substitution models (such as Q.plant or Q.insect). Is this likely to bias my data? Should I re-run the analysis and remove the substitution models that are meant for different groups?
Thank you for your time!

That is a complicated question @Laura and is certainly worrying. I do know insect mutation matrices and compositional bias but cannot think why that would be confused for bacteria. What is the %GC of the bacteria of interest? Secondly, what is the data, I presume it is nucleotide data.

Hello @M__ ! It's protein data, I forgot to mention. The GC is roughly 50-55%, with a few outliers (e.g. my outgroup has ~60% GC content).

Thanks only the protein bit matters ... and yes it is likely to make a different. Could you tell me whether it is using any other protein mutation matrix e.g. DayHoff, JTT or is Q.insect the matrix? It is fairly easy to generate a mutation matrix, rather than use an established one.

Thank you @M__    . The program is running the data on a lot of amino-acid exchange rate matrices (including JTT, Dayhoff, and a lot of variants of each accounting for rate heterogeneity across sites http://www.iqtree.org/doc/Substitution-Models). It seems to pick the best model on something called Bayesian Information Criterion? (it also lists the best model according to Akaike Information Criterion and Corrected Akaike Information Criterion but it doesn't look like they influence later steps as far as I can see). Overall (with the variants) it goes through about ~700 models each time.

... is there any chance the bacteria are insect symbionts @Laura?

IQTree is a fashionable method for speeding up bootstrapping and used on ultra large datasets. What IQTree might be doing is rigidly fixing a model and parameterisation is not estimate via the ML. This might be how it achieves its performance.
In context, IQTree can use an establised mixed model, which is one approach to protein phylogeny. As a personal note, I don't favour it and come back to the reason later. It MIGHT also generate its own mixed model. There are not 700 established (published) amino acid matrices plus established mixed models, combined there are no more than 50.
There are around 30 protein mutation matrices, Dayhoff models (I honestly don't who Dayhoff is), one for insects which might be for mtDNA ... mtArt, one for mtDNA, mtREV I think and so on. Please check carefully the number of models it is evaluating.
Its easy to see how 4 or 5 combinations (or possibly permutations) of 30 become 700 models if it is generating its own mixed model.
Alternatively, what it does (mixed model or not) BIC and AIC do is given a tree estimate (prior) what's model is the best probabilistic fit. In your case Q.insect is cropping up, whether this is mixed model or not ... dunno.
What I don't know is how it calculates the mixed model, is it creating a new singular matrix by averaging (which is probably closer to the algorithm), or is it sequentially passing the data through first one model, then another etc ... Personally, I avoid mixed models for precisely this reason, I dunno what the actual matrix is, but accept it can provide strongest AIC/BIC.
It suggests either:

BIC/AIC results in an outcome that isn't biologically realistic .... OR
there has been some weird parallel/convergent evolution - this happens in parasites for certain

If you are using a single matrix not mixed model and Q.insect is cropping up - I would be concerned. If it is part of a mixed model, I wouldn't worry.
Is it checking 30 models (plus published mixed models) or 700 matrices??
If the final point was correct (biologically real result), you check it via blast - simply exclude bacteria and see if you get a plant of insect sequence as the top hit.
Ultimately, the uncertainty is the reason I avoid mixed models.

To answer the question, I'm not aware of a bacterial amino acid matrix. HOWEVER, you can select a mixed model. It might just be screening around 50 models via AIC/BIC and within that Q.insect is the most probable.
Summary
The alternative is to look for support using a general amino acid matrix (for example JTT, LG or JTTDCMut).
Personally, I would be very cautious about proceeding with a matrix generated by insect genomes outside a mixed model for bacteria and in the absence of any biological support. FOR EXAMPLE, if these were insect symbionts - that makes sense.
Which ever approach is used a Bayesian approach is often required.


text:
How can I improve or otherwise investigate an unreliable genome tree?
Summary My genome tree doesn't agree with my gene trees and I get the feeling that my genome tree might be wrong, possibly due to long branch attraction, but I don't know how to check/fix it.
Background
I have a set of genes of interest recovered from genomes and metagenomes. I also have some bins built from the metagenomes of interest.
Previous work I've built multiple gene trees, which tend to agree with each other, as well as with other phylogenetic methods (16S rRNA, initial genome trees where available, conserved orientation of the genes of interest). So at this point I though I had a pretty good idea of how my genes / species of interest are linked to the species tree.
Concern
I've assembled the metagenomes and recovered bins. I've refined bins of interest with anvio, based on sequence composition / differential coverage, and kept only the bins which appeared quite high quality (>70% completion*, <10 - often <5% - redundancy). Thus contamination in the metagenome assemblies (MAGs) is unlikely. The contamination would have to happen multiple times in the same way to different samples from different environments, at the same coverage, to get the same effect.

with the exception of a bin of particular interest which I accepted with 62%

I have de-replicated these bins, added some known genomes, and tried to build a species tree with Orthofinder (with the proteomes). There is a problem occurs with a clade of 3 bins I will call 'X'.

Trees Based on the gene-of-interest phylogenies, X looks like probably a sister-group to the rest of known genomes from order 'A'. Which is very interesting.
With the 16S rRNA, it's not quite as clear (there are no 16S rRNA sequences in the bins); however, some 16S rRNA amplified from one of the original samples fit the bill as 'sister-clade to basically everything else in group A'). The first (100 or so) BLAST hits for this 16S sequence are all from group 'A', still.
Tree conflict However, in the proteome tree, X as well as Y (Y being a sequenced genome that I suspected was related to X based on 16S data) end up in order 'B' instead of 'A'. GTDBtk also labels X as belonging to clade 'B', not 'A'.
Now,I know that genes and genomes don't always have the same evolutionary trajectories. But multiple sets of data lead me to a different conclusion than the proteome tree, so I have to at least question whether the proteome tree is correct; what's likely to be the problem (I suspect long-branch attraction), and, if it is long-branch attraction (or anything else) how can I uncover the correct signal?
Thank you for your time!

Thanks @Laura ... I ain't sure LBR is occurring but would need to think. What is the tree building algorithm and parameters?

@M__ Thank you for having a look! I used Orthofinder (https://github.com/davidemms/OrthoFinder) which appears to use an algorithm called STAG (https://www.biorxiv.org/content/10.1101/267914v1) for the species tree building. I'm not sure on the details. The taxonomy tool GTDB-Tk (which appears to use HMMER /pplacer) is consistent with this proteome tree (but not with everything else).

Where is taxa 'Y' in the far left tree? I can see 'X' ... but no 'Y'

The difficulty in answering the question is:

The absence of "Y" in the "Gene of interest" tree
The sampling bias in the number of taxa mostly between the "genome tree" and the "gene of interest" and "16S" tree
The "genome" tree is unrooted so its difficult to compare it with the "16S" tree, thus its possible that if "16S" was unrooted it falls within monophyly "B"
The trees don't have any branch or statistical support
We don't really know the tree building method (I haven't looked it up)
The bacteria might have significant recombination (we don't know the bacteria).

The answer is there is a lot of room for different tree artefacts here, however it could be a biologically real result. I agree "the tree" is simply a tool which uses a series of assumptions and obtaining confidence in a tree requires a lot of analysis. Vast amounts of computational time are dedicated to obtaining a robust tree.
The principle step is to compare like for like:

same number of taxa and
same rooting position (or unrooted).
second step is to obtain statistical robustness for the observation.

There is a formal process for identifying conflicting trees and it needs code - I have written stuff that does this and maybe I should release it. Essentially, you would data chunk gene by gene to produce a genome of gene trees and use e.g. ETE3 or similar to plot the tree conflicts. It ain't rocket science. Thus the amount of user time is trivial. The idea is that once the "problem trees" are identified, e.g. rapid evolution.
The one observation I would say is that the "Y" branch looks a little bit weird and this could indicate part of the data is has resulted from different evolutionary rate and that is the basis for long-branch attraction.

Scenario 1 It could be that the 16S and "genome" tree are part of the same tree (monophyly B), they are just rooted differently. Thus the conflict is just a small number of genes and this could lead to ...
Scenario 2 Keep in mind bacteria do undergo 'horizontal transmission' and if 'gene of interest' is under strong selective pressure recombination is entirely possible.

Hello. My mistake - I get these results with the rooted trees. HGT would be possible, but hasn't been shown for that particular set of genes before, so I would have to provide fairly strong evidence. I also think it would be highly unusual (assuming that organism X really belonged to B) that when blasting the 16S of an organism, the results I'd get would mostly be from another clade. Thank you very much for clarifying this and telling me about ETE3. The program I used for the genome tree has the gene trees stored, so I can compare them : )

Great @Laura glad its useful and its worked out. ETE3 toolkit BTW, bit different from ETE3 python module is stuck in Python 3.6. The all-singing all-dancing ETE4 was supposed to be released already - it got full funding, but maybe never.

It seems like you've got the right approach, in trying out a few different genes and looking for consistency. As you say, different genes can have different selective pressures, which influences how much they change in relation to the whole genome background.
Certainly don't put too much trust in a single gene, and especially not 16S: it has a low mutation rate, can jump around and be repeated in the genome, and is an occasional target for horizontal gene transfer.
The answer may be that there isn't any single correct phylogeny. A primary assumption of phylogenetics is that evolutionary progression has a tree-like structure (with roots and tips), but that model is broken by horizontal gene transfer and plasmid sharing. It may be the case that you need to break the genome up into fragments, each with its own inherited history.


text:
What does the 'suppressed' folder mean on the NCBI ftp genomes website?
I wish to understand the term suppressed in folder nomenclature on the NCBI ftp genomes website?
Example: here ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Oscillatoriales_cyanobacterium_C42_A2020_001/all_assembly_versions/
Normally, I'd see folders with names like representative or latest assembly, which are pretty straightforward. But what does suppressed refer to in this case?
Thank you for your time!

Here you have the explanation:

Suppressed genomes: Several hundred RefSeq bacterial genomes that did
not pass assembly or annotation quality validation have been
suppressed(RefSeq genome quality criteria). An FTP report file listing
the most recent suppressions will be provided soon in the RefSeq FTP
site.


text:
Are there any resources for calculating / representing evolutionary relationships within gene clusters?
I am working with a gene cluster found in several species. The genes involved are pretty much the same, with the odd insertion, but they can be often found in different orientations. So, I can't just do a multiple alignment of the whole cluster. I've done phylogenies of individual genes and I've looked at cluster architecture (indels, rearrangements) and things make sense (eg. limited HGT), but..
...Is there a simpler way to do this? A way to calculate / represent a phylogeny of the whole cluster? Like...some sort of web, maybe? Most of the software that does take rearrangements into consideration is meant for whole genomes, and is unreliable when it comes to fragments about 20-30kb long (like mine).
Thank you for your help!

Can you give us some example data? It's kind of hard to understand what you're looking at. Are you trying to build a phylogeny based on the entire cluster or have you separated each species' cluster into a multifasta file with each gene separate? If you haven't, wouldn't that be a reasonable approach?

If there'd be a way to build a phylogeny based on the entire cluster (that can deal with indels/inversions/rearrangements), I'd love to know it.
I guess I could separate each species' cluster into a multifasta file, but I'm not sure what I could do after that? Could you please explain? It's unfortunate, because I'd lose positional information (eg. let's say species X is closely related to species Y and less so to species Z. In X and Y, the genes are in the order A-B-C-D, in Z they're in the order A-D-C-B), but if it's the only way...

Well, I don't really think the order will make much of a change in the phylogeny. I'd just use the multifasta approach to build the tree and then maybe annotate it with the directionality.

Its easy ... I strongly suspect these are bacteria ...please confirm. The likely application is resistance cassettes again please confirm. Finally this is likely a "genetic island" rather than a single gene ... again please confirm.

@terdon - Thank you. I've never built trees with multiple genes at a time, before. Could you recommend some software for it, please? Until now I've mostly used Seaview and the CIPRES webserver.

@Michael G. - Bacteria, yes. (Silly me, forgot to specify); resistance cassettes, no. And, yes I'm talking of a group of genes rather than a single gene.

How many taxa? 1000s or <100 that is kinda important. I get the idea the fragment is 20-30 kb

@MichaelG. Around 20 taxa.

Oh that's easy. I'll post a full response later/tomorrow if no one else does.

Firstly you wanted a webserver to build a tree, ok, I would use PhyML server. The good thing is that with only 22 taxa the calculations are easily doable.
What I'm assuming is you can map your rearrangements to produce a single homologous gene alignment for all taxa. So you appear to infer you can break your 20-30kb contig into individual aligments of homologous genes which has incorporated the rearrangement and enables easy production of an alignment file. This would leave you with around 20-30 gene alignment files. I hope that is correct, because that is how you get around the rearrangement problem and enables the data to be reassembled for the final tree (below).
Doing the tree is simple, because you can exploit the assumptions used in tree building regarding independence between mutations at each alignment position. The theoretical proof, the data is sound, is complicated. In summary,

build a consensus phylogeny (one tree representative of lots of trees) using two different approaches
combine the data into a single alignment, if the above looks good, or excluding the alignments that are affected recombination.

The theoretical issues with what you want to do are complicated because you need to distinguish bifurcating (clonal) homologues from recombination that could have occured in the reassortment events. The following is a  purists approach, you might not want to do this stuff, but it gives you some idea.
Just before I do that ... I just want to check the alignments assumption is correct? and which bacteria?


text:
Variation in 16S rRNA between assemblers - how do I know which is real?
I have a low-diversity metagenome (~11 bins > 80% completion). Out of the bins, 3 are of interest to me. None of the lower-completion bins that can be identified are from the group of interest. So let's call the bins of interest 1, 2 and 3.
The raw reads had been assembled separately with MEGAHIT and SPAdes (meta). Out of them, only the SPAdes assembly was binned / refined. There are 3 16S rRNA in the MEGAHIT assembly (x, y, z) and 3 16S rRNA in the SPAdes assembly (a, b, c).
Bin 1 is associated with sequence 'a', which is identical to 'x'.
Bin 2 is linked with SPAdes sequence 'b', which has no MEGAHIT equivalent. Bin 3 is not associated with any 16S rRNA (even with MarkerMAG).
Sequences 'c' and 'z' are very similar. Sequence 'y' is an outlier.
There are small hints from phylogenetics that bin 3 could be associated with sequences c/z. But then, could sequence 'y' just be a misassembly? Or is it that somehow my final bins are missing an entire genome of interest? If the 16S rRNA for 'y' is real, I would expect a genome divergent enough from the others that it wouldn't get binned together with them. This is what I am trying to understand. Would you happen to have any advice / suggestions on what I could to to test this?
Initial approach: assessing scaffold quality with the BBMap stats.sh . It's a bit difficult to tell, since SPAdes has about 10x more contigs. Worse N50/N90 but more long contigs and a better 'general assembly score'.
Additional thoughts: The new bins seem too divergent from known genomes to use a reference for misassembly estimation.
Thank you for your time!

This is easy @Laura - I'm a bit surprised by the results Spades meta didn't have a good reputation, but the criteria assessed was not what you are looking for.
Anyway take all outputs and place them in large multi-species alignment,

e.g. for each output perform a Blast for maximum genetic diversity for each sequence and sample across the large Blast  output. Take note if there are multiple partial hits against diverse bacteria (online Blast has a diagrammatical summary).
You have a good understanding of your taxonomic group and intersperse that with other stuff, i.e. fill in the diversity.
Assemble it, e.g. muscle5 superfast option.
Perform a split decomposition phylogeny or equivalent which identifies phylogenetic networks.

SplitsTree is the historic package for this and it turns out that split decomposition was not the best, I think "neighornet" was superior. What this is does is identify mosaic gene signal in your data. It provides a visual representation. Thus if your data is a artifactual mix of two or more taxa this will cause a network effect between the two parents and your newly mis-assembled sequence will sit between the two in a two "diamond same". This is the antonym of a bifurcating phylogeny.
The final step is the assembly that produces sequences that minimises phylogenetic networks is the true sequence. If you've not encountered split-decomposition/networks before a verbal explanation might seem confusing but when you run through a few examples and try it on your data you'll get the grasp.

Much appreciated! Good idea. Tested. Couldn't see obvious signs of a chimera, but the alignment showed that most of the short 'outlier' sequence's distinguishing features were not present in any other (even distantly related), including a lot of weird Gs. Chimeric with something else maybe? Curious thing is, the depth of that MEGAHIT contig seemed to be pretty good (~120x). I had thought I'd done my homework by quality-controlling the reads with BBDuk, but it's a bit unsettling to see how software can be so certain and so wrong at the same time. Should look more into missassembly detection...

It tentatively sounds like you might have found a new species. Just keep in mind the diversity of bacteria in the analysis affects the ability to detect a possible chimera, i.e. if the "parents" are not present it could throw a false negative -------- Weird "Gs" high AG are a known feature of bacteria: particularly if its genetically distinct. Anyway if Splitstree if done correctly it's a sensitive method for chimeras: but you might have to shift the sample size - remove distant sequences which can cause a collapse in networks.  Megahit is respected BTW.

It's more like...The positions in questions are 100% conserved within related bacterial orders, and here's this short sequence with ~9 of these otherwise conserved positions being different. I guess I can't prove it's not real (part of my issue...), but it doesn't look...correct?

What you are saying sounds perfectly reasonable - keep in mind I can't see the data nor analysis. Proof ... it's the first stage requiring downstream analysis, so either leave it as a "putative result" or go into the complexity of metagenome assembly by assembling other single copy gene targets like atpD.


text:
Is there a graphical/interactive 16S rRNA clustering method?
I've been doing phylogenetics with large (hundreds) 16S rRNA sequences lately.
Usually I'm focusing on one order, and using a combination of trees and sequence similarity to assess stuff like 'is this a genus', 'is this a separate genus', 'is this sequence actually in the order of interest'. It feels like a lot of what I do is repetitive and could be automated more. In particular, it feels like there should be a simpler way for having an assemble view at sequence similarity rather that calculating a Percentage Identity Matrix and then putting it in Excel with conditional formatting. I'm imagining something like this clustering diagram for species (https://www.researchgate.net/publication/360738216/figure/fig4/AS:1182363421806594@1658908815491/The-t-SNE-clustering-of-bacterial-genomes-in-the-code-RA-space-at-different-taxonomic.png), but for sequences instead (using sequence similarity values, at least). Still, I can't seem to find something like that for 16S rRNA that's as visually compelling (or even as visually helpful as my Excel sheet). There are a number of programs for clustering OTUs, but the ones I've encountered don't seem to have a graphical / interactive interface.
Does something like this exist? Thank you for your time!

My recommendation is MegaX. It's the closest thing in phylogenetics to an all encompassing GUI.

This has a nice visual alignment editor which highlights conserved, variable and parsimonious sites.
It will produce distance matrices based on any number of criteria, e.g. % homology, Jukes-Cantor, Kimura-2 parameter etc ...
Good for bootstrapping distance-calculations as well as can perform ML based tree building methods with full parameterisation.
Nice tree visualisation, which can be edited, switched to viewing bootstrapping (if performed), or to branch genetic distances (again under any distance criteria, e.g. %age homology etc ..). It does loads of other stuff, which in this case are not needed.

You can export the tree and import into FigTree to produce collapsed clades represented by coloured triangles. In summary for easy GUI stuff MegaX, then FigTree.
There are no PCA/MDA based analysis for 16S sequence data because eiganvalue-based calculations would override lots of theory involved in point mutation modelling and the resulting trees. I could replicated the PCA output in the question with sequence data, its not hard to do - but it's not publishable.

Thank you. I'm not sure I understand the last part. (I'm not really familiar with PCA/MDA/eigenvalue) Is it difficult for some reason to represent sequences as blobs, with distances between them based on sequence similarity values? (and potentially different colors based on tresholds)? I don't think that sequence similarity alone is a good idea, but it has provided helpful supporting data to trees in the past (e.g. 'yes, these genera are likely closely related even though support values are abysmal' or 'no, there is no support for these genera being closely related, this is probably LBA')

@Laura, the 'difficulty' is the diagrams you showed were a PCA analysis (axis were 1+2 components). The problem is not the practicality but the theoretical model of what a 'blob diagram' represents. The theory you are using is modelling point mutation, Jukes-Cantor correction and all that. You can't replace Jukes-Cantor with PCA - completely different.  Also you can't do PCA on a distance matrix because there's no multivariates. Sure, I could code a 'blob diagram' but the theory used would be controversial - thats the problem.

Converting a distance matrix to a 'blob' diagram - it's better to use clustering analysis, i.e. make a tree. Clustering algorithms (you normally use Neighbor-joining but it's one amongst many) are the best approach to analysis of a distance matrix. LBA + support ... if could put that as another question might be able to help.

We can chat briefly here https://chat.stackexchange.com/rooms/58858/bioinformatics but I will simply explain that the model used isn't applicable to blob diagrams

Thank you. Sorry. I lack the vocabulary. I guess I wasn't really aiming for a PCA, just a graph where I could show that different points (in my case, sequences) are similar by how close they are together. Sort of, sequence A and B and C would be close together due to high similarity values (~97%), and D would be a bit off to the side (but closer to C than A or B). I'm not sure if this makes sense. I'm also confused as to how sequence similarity values necessarily imply an evolutionary model or evolutionary distance (blastn has a scoring matrix of +1match/-3gap, but does that imply J-C? @M__♦


text:
How do I compare transcriptomes of the same species across different developmental stages?
I want to identify the most highly expressed genes across several developmental stages of a given insect, in hopes of identifying a protein that can be used to design an insecticide.
I have a few doubts :
Objective

Several studies have been performed on individual developmental stages and their transcriptome has been published. Should I compare each set of data against the other manually, to find genes that are expressed universally and then target the most highly expressed ones or is there any software / computational method that can be used instead?

Question
Is it right to compare transcriptomes to find commonly expressed genes?
I did try to look at the materials and methods section of a few papers on comparative transcriptomics, but it only says that data was compared without stating the method. One paper cited VennBLAST, but I wasn't able to find the tool online.

Could you clarify what do you mean by individual developmental stages? They looked at the same individual in several stages or they look in several samples different stages of development. What I think you want to do is find differentially expressed genes, search for those terms and let us know if you need something else

I mean that each paper concentrates on a single developmental stage. They analyze several samples for different stages of development.
Also, no, I do not want to find differentially expressed genes.From what I understand, differentially expressed genes are those that are upregulated only in a particular cell type (in my case stage). I want to find out if a common set of genes are being upregulated throughout the life cycle of the insect, ie. they are common to all the life stages. I thought I can do this by analyzing the transcriptomes of each stage. I hope this is clear!

You can find the differentially expressed genes that are downregulated too, and it is not specific to a cell type that depends on your data. But it is a comparison between two stages, so you can compare all the initial stage with all the following stages. Then you find the common genes that are up-regulated in all the cases of interest

Oh! I didn't know that differentially downregulated genes can be found too! Thanks for your answer :D

This was a comment!! You seem to need more guidance, ask your advisor or a senior people how to find the most expressed genes and from there how to find it those that are so in several stages.

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


text:
Meaning of category in MutSigCV
I am trying to understand the working of MutSigCV. But I am not sure about some terms used in MAF file format which is category.
According to MutSigCV, every mutation can be divided into the following categories:

CpG transitions
CpG transversions
C: G transitions
C: G transversions
A: T transitions
A: T transversions
null+indel mutations

I am very much clear about some of this categories like C: G transitions and A: T transitions but for rest, I am unable to make any solid conceptual background. Can anyone suggest me the proper definition of these categories of mutations with proper examples?
Thanks.

Definition
MutSig stands for "Mutation Significance" and is defined here

Did you try asking on the help forum? Or your problem is what is each mutation?

@Llopis: Not yet. I even didn't know about that. I know only GATK forum and this website. I have googled these terms extensively but didn't get even proper definitions. That's why I started the bounty on this question.

CpG transversions are the same as C: G transversions, except that they are followed by a G which means they are potential targets for DNA methylation. Methylcytosine can undergo deamination to Thymine, which makes them mutation hotspots

indel = insertion/deletion. I have no idea what null means in this context

@heathobrien Those comments are an answer! With a link to a source you could get the bounty


text:
Post ipyrad filtering of SNP loci from GBS, based on SNP quality
I want to use SNPs produced by Ipyrad, which is a python script for RADseq, using maize genome data via RAxML to examine the monophyly of a highly-variable focal species and its phylogenetic relationships with about 40 congeners.
Essentially phylogenomics using maximum likelihood trees.
An initial analysis suggested that samples from two outgroup species might share few SNPs with other samples, so the threshold for shared loci was set to two samples.  As a result, the .u.snps.phy alignment has over 600,000 SNPs.
For subsequent use in RAxML, I would like to implement two filters to construct alignments with subsets of SNPs:

SNPs that occur in at least one outgroup sample and in some (potentially maximal) number of other samples,
SNPs constituting a random sample from an arbitrary subset of samples, i.e. potentially excluding the outgroup samples or other samples,
merge the alignments into one

How can I do this efficiently, maximizing use of existing functionality in R or Python? Either code or a (likely successful) strategy would be helpful. I have access to all output formats from ipyrad as potential starting points.

Sorrey but I was lost at the first sentence, could you add a link to what ipyrad, GBS and RAxML is? Also, could you explain what have you tried ? (It is better to show the effort done before asking random people in internet) You seem to have already a plan, so you might have already coded something...

My bad.  Links are now available. The initial analysis wasn't appropriate and I didn't code it. I have no code yet, or I would have posted it. I do have some idea about what I want the code to accomplish. Since I haven't post-processed an alignment before, I aimed my question towards a non-random someone with knowledge of ways to post-process alignments (and specifically from ipyrad).  I naively tried a tibble (in R), with SNPs in columns, but it hung--must not be an appropriate approach with this much data. If I could work on rows and cols, it would be easy, but another approach must exist.

I don't know the answer but have a look at Bioconductor packages (R) and Biopython (python).

Its easy, just Biopython's AlignIO, all manipulations can be done therein. RADseq is restriction digest genome sequencing, bits of random contigs. I think the question is an example of maxium obstrufication, its quite simple really..


text:
PAM50 gene expression classification
I'm looking at running the PAM50 classifier on RNA-Seq from 138 breast cancer samples. However, the R package (genefu) that's useful for this does not have a particularly helpful vignette when it comes to processing the RNAseq counts because it uses pre-processed publicly-available RNAseq data.

Starting from counts, is there a way to pre-process the data for running intrinsic.cluster.predict in R/genefu using the built-in signatures? (I probably won't have enough power to split the set into a training/test set).
The best choices for the parameters in intrinsic.cluster.predict would also be very helpful.

For differential expression (I have pre and post surgery pairs), I have used conditional quantile normalisation to get normalised log2(RPKM) values for running through R/EdgeR, so I also have these available if that is the best way to go.

EDIT: looking at the genefu vignette, this is all based on microarray data (which is essentially log2(expression)+offset), so will adding an offset to the log2(RPKM) values improve prediction for PAM50?


text:
Repeat analysis on eukaryotic assemblies
I have a hundred insect genomes and I'm looking for repeated regions along these assemblies.
First of all, I thought of using RepeatModeler, only the most familiar among you know that it is not at all optimized software that requires a lot of resources and calculation time, which is therefore not feasible in my case.
One possible solution Repbase is not free.
I thought for example to use a blastn approach by using in query the assemblages and in database of nucleotide databases containing ETs, which would be an approach maybe less sensitive but much faster! Only I can't find a database in ET representative of eukaryotes
I have two questions:

Is there an alternative algorithm that is fast enough to process a large number of insect genomes?
Is there a database of representative ETs to enable a blastn style solution?

Thanks a lot.

when you say "ET" do you mean "TE" ("transposable element")? If not, what does that abbreviation stand for?

Yes Transposable element sorry