text:
how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)
i have a question when i want to extract cds sequence using gene id. but cds file is not just start with >gene is, it has many other annotation. the only same  is star with gene:
cds fasta:
>Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding
ATCGGCAT
>Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA

gene_list
Zm1d035916
Zm1d048

how to get result like
>Zm1d035916
ATCGGCAT
>Zm1d048
ATGCGGCA

If you can use awk, you could use this simple script:
awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

You can also put the code inside a script file:
extract_genes.awk
#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

And call it like this:
awk -f extract_genes.awk gene_list cds.fasta

hi, thank you for your answer. Where to find my output file? After I run it, I don't know where to find the output file when finished. Thank you!

Happy to help. The output is written to stdout. You can forward it into a file like this: awk -f extract_genes.awk gene_list cds.fasta > output.fasta

Thank you. However, I don't know why the output.fasta file is empty. I use command like this: awk -f gene_cds_extract.awk gene_list_red_green_cpsf100.txt all_species_longset_cds.fa > example_cds.fa

That's odd. Are the input files formatted the same as in your question? And are there overlaps between the gene names and the fasta headers?

Thank you for your quick reply. input cds.fast format like this: >PNS24249 cds supercontig:Brachypodium_distachyon_v3.0:Bd1_centromere_containing_Bradi1g41430:24784:29531:-1 gene:BRADI_41430s00200v3 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGCGGCCGCGCGGCAGAGGAGACGACGGGGACGACGACGAAGACGAAGAGCCCTCGCGT
CAACGGAGCGCGCCGGAAGCCGACCCTGATCCCGAGCCCGCGCCAGCGCAGCAGCCCACA
CGAGCGCCCCTCAGCAGCCTAGTCGTCAGGCCGCCGCCGCCGCAGGACAACGGCGCCTCC     input list.txt: BRADI_41430s00200v3 like this. no gene overlap, I extract longest cds in cds.fa file

I just made a test with your input and the script above and it outputs the sequence. Which operating system do you use, and which awk version (I use Lubuntu 18.04 and GNU Awk 4.1.4)?

One thing you could try is to use a simpler (less specific) regular expression. If you replace re = "\\sgene:"g"\\s" with re = "gene:"g does it output anything?

I use mac terminal. I don't know it already have awk or not. Sorry, It's my first time to use awk in terminal.

Thank you, I try it with simpler regular expression, and it got results. Great!!! Thank you very much.

Aha! Great that it worked. Maybe awk on Mac understands special characters in regular expressions (in this case "\s" for a space) differently than awk on Linux.

Anyway, it saves a lot time to extract sequence, not just find, copy and paste, you make it more interesting. Thank you!

If you read it in python, they would take care of that:
from Bio import SeqIO
records = list(SeqIO.parse("test.fa","fasta"))

records

[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm002', name='Zm002', description='Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm001', name='Zm001', description='Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding', dbxrefs=[])]

Get the gene and sequence.. maybe there's a easier to do this in biopython, this is like my 10th time using it:
def format(rec):
    des = rec.description.split()
    gene = [entry[5:] for entry in des if "gene:" in entry][0]
    return(">"+gene+"\n"+str(rec.seq))

[format(rec) for rec in records]

['>Zm1d035916\nATCGGCAT', '>Zm1d048\nATGCGGCA']

Write it out:
f=open('example.fasta','w')
f.writelines("\n".join([format(rec) for rec in records]))
f.close()

list(SeqIO.parse("example.fasta","fasta"))
[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm1d035916', name='Zm1d035916', description='Zm1d035916', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm1d048', name='Zm1d048', description='Zm1d048', dbxrefs=[])]

This is a very neat comprehension

lol OP chose the mighty awk =p . Thank you for the kind comment

The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records]))

Thanks! I try it and it is also work.


text:
How to map gene names to ensembl Ids in R
I am trying to convert a list of gene names into Ensembl IDs but I'm having a bit of trouble.
I've tried using the EnsDb.Hsapiens.v86 package but I always get LRG-ID values returned as well, even when using the GeneIdFilter "startsWith" thing.
My code looks like this:
targets <- as.character(targetCSV[, 2])

targetEnsemble <- ensembldb::select(
      EnsDb.Hsapiens.v86,
      filter = GeneIdFilter("ENS", "startsWith"),
      keys = targets,
      keytype = "SYMBOL",
      columns = c("SYMBOL","GENEID"))

And for example I get back:
    SYMBOL          GENEID
1     CYR61 ENSG00000142871
2      CTGF ENSG00000118523
3    AMOTL2 ENSG00000114019
4    ANKRD1 ENSG00000148677
5    ANKRD1         LRG_379
6    IGFBP3 ENSG00000146674
7        F3 ENSG00000117525

I also tried using org.Hs.eg.db package but get a different problem where some of my genes return NA.
My code:
  mapIds(x = org.Hs.eg.db,
     keys = targets,
     "ENSEMBL",
     "SYMBOL",
     fuzzy = TRUE,
     multiVals = first)

And this works on some but not all the genes for some reason.
    $CYR61
[1] NA

$CTGF
[1] NA

$AMOTL2
[1] "ENSG00000114019"

$ANKRD1
[1] "ENSG00000148677"

$IGFBP3
[1] "ENSG00000146674"

$F3
[1] "ENSG00000117525"

If anyone can solve either of these problems that would be great, thanks in advance!

Have you considered just downloading a list of gene IDs and the associated gene names from Ensembl Biomart?

Yeah but I wanted to be able to change out the genesets quickly and this just added another annoying step

There's an R package to allow fetching those programmatically from biomart.

source: Stack Exchange

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

sort: votes


text:
Ranking metric in GSEA
I am using fgsea in r to calculate and plot a bunch of GSEA graphs.
My question is whether it is acceptable to use the Wald statistic from DeSeq2 to rank the gene list?
I have seen in the GSEA application that the signal to noise ratio is used but is the Wald statistic sufficient?
If not then is there an easyish way to calculate the signal to noise ratio in R given a normalised count matrix?
Thanks in advance.

If using DESeq2 with GSEA, I'd recommend ranking by shrunk log2FC values.
It'd also be worth considering ranking positive and negative associations separately, because the standard GSEA algorithm doesn't cross at the zero point when associations change from positive to negative.
You shouldn't be using p-values to rank anything. A p-value gives information about the reliability of a result, but not its importance. See point #5 of this article about p-values:
https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i31

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


text:
ERROR: (subscript) logical subscript too long
I am trying to remove a row of a dataframe in R however because the dataframe is so large when I use:
rawCounts <- rawCounts[rawCounts != geneID, ]

I get the error:
ERROR in rawCounts[rawCounts != geneID, ]: (subscript) logical subscript too long

Is there another way of removing a specific row in a large dataframe that is not so memory intensive?
Thanks in advance.
I have quite a bit of memory, memory.limit() returns:
> memory.limit()
[1] 16318

I am answering my own question here.
I realised I should have been using
rawCounts <- rawCounts[rownames(rawCounts) != geneID, ]

Problem solved.


text:
how could i find uniprot ac in blast?
I have made a search in p-blast.
I would like to search for the protein results in uniprot. how can I do it directly?
I can go to GENBANK but i would like to go to UNIPROT.
Thank you!!!

Ah ... yes thats tricky. The only thing I can suggest is using Biopython for Genbank and Uniprot because Biopython facilitates both.

thank you a lot for your answer!!

UniProt also provides a BLAST service (https://www.uniprot.org/blast), or you could use the Expasy BLAST too, at https://web.expasy.org/blast.
If you have performed your similarity search at NCBI or elsewhere, you can map the identifiers from that result to UniProt by using the UniProt IDmapping service at https://www.uniprot.org/uploadlists. This can also be used programmatically, cf https://www.uniprot.org/help/api_idmapping

This is all very useful information... just want to add that NCBI Protein BLAST also has "UniProtKB/Swiss-Prot" as a Database option.


text:
can I download DESeq2 in R 3.6.3 in Linux MInt?
I am trying to download DESeq2 in R 3.6.3. Is it possible?
This is printed when I am trying: Warning in install.packages : package ‘DESeq2’ is not available (for R version 3.6.3)
Thank you!

You might have more luck using conda to look for older versions of the DESeq2 package that are compatible with R 3.6.3. If not, you can use conda to install compatible versions of R and DESeq2 in a new environment and use that environment for your analyses.

Bioconductor is a package management tool for R that allows you to install the appropriate package versions for your R version.
For R 3.6.3 you'll need Bioconductor 3.10:
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.10")

Then you install  DESeq2:
BiocManager::install("DESeq2")

these errors are generated when I am doing what PPK suggested:
Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3
(2020-02-29)
Installation path not writeable, unable to update packages: nlme,
spatial, survival

BiocManager::install("DESeq2")
Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3
(2020-02-29)
Installing package(s) 'DESeq2'
Warning: dependency ‘XML’ is not available
also installing the dependencies ‘annotate’, ‘genefilter’, ‘geneplotter’

trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/annotate_1.64.0.tar.gz'
Content type 'application/x-gzip' length 1873132 bytes (1.8 MB)
downloaded 1.8 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/genefilter_1.68.0.tar.gz'
Content type 'application/x-gzip' length 1421076 bytes (1.4 MB)
downloaded 1.4 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/geneplotter_1.64.0.tar.gz'
Content type 'application/x-gzip' length 1431450 bytes (1.4 MB)
downloaded 1.4 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/DESeq2_1.26.0.tar.gz'
Content type 'application/x-gzip' length 2060061 bytes (2.0 MB)
downloaded 2.0 MB
ERROR: dependency ‘XML’ is not available for package ‘annotate’

removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/annotate’
ERROR: dependency ‘annotate’ is not available for package ‘genefilter’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/genefilter’
ERROR: dependency ‘annotate’ is not available for package ‘geneplotter’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/geneplotter’
ERROR: dependencies ‘genefilter’, ‘geneplotter’ are not available for package ‘DESeq2’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/DESeq2’

The downloaded source packages are in
‘/tmp/RtmphLXrZb/downloaded_packages’
Installation path not writeable, unable to update packages: nlme,
spatial, survival
Warning messages:
1: In install.packages(...) :
installation of package ‘annotate’ had non-zero exit status
2: In install.packages(...) :
installation of package ‘genefilter’ had non-zero exit status
3: In install.packages(...) :
installation of package ‘geneplotter’ had non-zero exit status
4: In install.packages(...) :
installation of package ‘DESeq2’ had non-zero exit status
I have tried to instal XML and the other dependencies without success..


text:
why in RNA seq don't we only use reference transcriptome?
I would like to ask why in RNA seq analysis (alignment step) we use sometimes reference genome instead of reference transcriptome? thank you!

Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.
But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.

thank you! Could I use STAR or a similar program (not Kallisto or Salmon etc) with reference transcriptome?

You can make a STAR reference out of a transcriptome, but why would you do that?  It wasn't designed to be run like that.

if i am only interested in differential expression, could i use star with reference transcriptome to decrease computational effort?

As a rule, do not use programs for things they were not designed for unless you are absolutely sure you understand what you are doing and why.  If you want something faster than STAR, use Kallisto.

I actually think you could run STAR on a transcriptome. In their manuscript, you can run alignment using a GTF file.

I didn't say you couldn't, just that there's no good reason to.  The point of a GTF file is to mark features, especially exons,  in a large genome, so there's little point in using one on a transcriptome.  And of course you can't align to a GTF alone.

I'm confused, the GTF file is the reference transcriptome no? And the FASTA file the reference genome (of course both are needed to map to transcriptome)

@Chris_Rands I wouldn't have spontaneously called a GTF file a transcriptome. For me, a transcriptome would rather be a set of sequences (the sequences of the potential transcripts). But maybe this is not the common usage.

@bli I agree, the transcriptome is maybe a set of mRNAs and not simply a tabular file format of course, but I was thinking in terms of how this information is stored

@Chris_Rands I never worked with a transcriptome, but I imagined that this was stored in the form of a fasta file containing all possible transcripts. And I imagined that mapping on a transcriptome was done via an index built on this fasta file. Is it actually more common to map on a genome and then infer transcripts based on information present in a GTF file?

@bli i think splice aware aligners need the full genome FASTA but others use transcripts FASTA


text:
why in 3' quantseq should we use reference genome and not transcriptome?
I am studying about RNA seq. I have found that in 3'prime rna seq (quantseq, tagseq) we should use a reference  genome  with    good    annotation  (plus known UTRs). My question is why?
Why could not we use reference transcriptome?
Thank you!

You asked a similar question before https://bioinformatics.stackexchange.com/questions/14093/why-in-rna-seq-dont-we-only-use-reference-transcriptome and were told to use the genome. If you so badly want to align to the transcriptome then do that but stop asking us for advice. Asking something when you are already convinced of an answer is a waste of time. Repeating the question when you don't like the answer is just plain rude.

It would be great if you can elaborate on where you "found" this information.
In the transcriptome annotation you have a defined number of transcript per genes, but there might be others that are present in your sample and not in the annotation. For example you have treated the cells with some extreme condition that changes the 3' processing machinery.
With 3' quantseq, you are left with the read fragment that is at the end of your transcript. If the ends of these transcripts are not part of the sequences in the annotated transcriptome, these will be unaligned.
These instances might be rare if your transcriptome annotation is good, but let's say if you find that a gene of interest which should be expressed is not, and you suspect it might be due to a different isoform used, you can always check it with the genome alignment (see where reads might be falling in the gene). With a transcriptome alignment, this will not be possible.

I'm not sure that's true:
https://www.nature.com/articles/s41598-019-55434-x
"QuantSeq. 3′ Sequencing combined with Salmon provides a fast, reliable approach for high throughput RNA expression analysis"

Yes I have seen this paper..thank you..but on another reliable source I have found the piece of information that I posted..


text:
two aligners combined results
I would like to ask if someone uses two different aligners to produce count matrices and then run the same alogrithm for DEG analysis, would it make sense to find the intersection of the DEGs in order to conclude to the "real" DEGS?

Thank you so much for your feedback!Yes that particular problem I am trying to solve..!

No, I am just trying to understand which algorithm to trust and which approach would give the most reliable results.. I am not trying to tweak results..

This sort of procedure is done in tool papers to benchmark results compared to different aligners/mappers. See the salmon or kallisto papers for examples. Using the intersect of DE results from different mappers will end up yielding somewhat conservative results, not unlike using different differential expression packages and taking the intersect of their outputs. In practice you could just take a more conservative p-value cut-off and save the hassle of using multiple tools, the results would end up being similar.

Yes, I have read in tutorials that the best approach would be to find the intersection of DEGs produced by different algorithms (deseq2, limmavoom, edger) and that s the reason why I wondered if I should the same with the aligner as well.. Thank you for your answer!


text:
how can i search in p-blast against UNIPROT(not only swissprot)
i need to search for one protein against UNIPROT in BLAST for a selected organism.
How can I do it? (SWSISSPROT+TREMBL)
THANKS!!

I don't know of any online tool enabling a search against TrEMBL. You probably will have to set up your own local BLAST environment, roughly following these steps:

download blast+ - but there might be an easier way to install the software, depending on your operating system (e.g. on Ubuntu you can just sudo apt install ncbi-blast+)
download target databases
build local database

I hope this helps.

thank you a lot for your answer!unfortunately my programming skills are not so strong!

UniProt provides a BLAST service (https://www.uniprot.org/blast) which you can run against several UniProt databases, or you could use the Expasy BLAST too, at https://web.expasy.org/blast.


text:
how to evaluate neural networks
I have made some architectures of neural networks in order to classify proteins in three categories. I have calculated true positive, true negative, false negative, false positive, Mathews, sensitivity, specificity, accuracy for each category. My question is what is the right way to evaluate the architecture. Except for the c-index, could I calculate a cumulative sensitivity, accuracy etc?or should I evaluate separately for each category? Thank you!

Wow, I would just go for accuarcy as the singular most important measure, but be careful of sampling bias.
If the classification training had an unbiased sample* and you have an accuracy of >90%, in general ANN that wouldn't be great, but in bioinformatcs thats really good. If your training comprises sampling bias that gets tricky and you either need to accomodate that in your accuarcy ROC etc ... or else retrain to eliminate the sample bias.
You can account for all of them of course by wrapping everything into a single analysis.
*, by this I mean that each of the 3 classifications comprises 33% in each training set. If this is badly skewed, this requires consideration and gets complicated. If it is precisely 1/3 each then 'accurary' is the key stat.

thank you!how can i calculate accuracy for all the three categories? is it correct to calculate the mean ? (sorry in advance for this question)

Thanks, firstly this looks very cool work, again... its cool.  Did you do 3 independent training sets, each with a single classification OR a single training with 3 classifications?

a single training with 3 classifications.thank you!

I get it, sorry. Yes and no. Yes: for your own internal purpose yes and a line or two in the report/paper. You don't want major skews between the classifications (these are different skews from the sampling skews). No: for the overall accuracy score, its one overall number (obviously one which isn't hiding a skew). Essentially the answer is its just one number, but you need to be confident its representative. I hope thats ok.

If you get 90% accuracy with a 33% sample for training of each classification, then its taken as read that the misclassification has reasonable distribution. Leading bioinformatics publications have present been as low as 80%.

so to sum up, instead of  making a 3*3 matrix to measure accuracy, I could calculate the mean (because now I have constructed matrices for each category). and if this is higher than 90% (but even less) then my model is good. am I correct? Thank you so much!!!

Basically yes in my opinion 90% would be fantastic and you are sorted, look at past papers (Bioinformatics journal, to verify the accuarcy against their work ...  Your confusion matrix (the formal name for a 3*3 matrix) helps diagnose your training, if your accuracy is low, scikit-learn has an automatic function for this. Personally I think 85%+ is good, but just my opinion.

thank you!in order to select one architecture than another how much should be the difference in accuracy? for example from 0.856 to 0.860 is it considered imortant difference?from 0.85 to 0.86?or it needs to be even bigger lets say from 0.85 to 0.90 or more?thanks!!

Either 0.856 or 0.86, sounds like what you really have, I think is fine and shows archecture convergence. Please refer to other comparable papers however for accuracy comparision. Obviously 90% is better. Things get very complicated indeed increasing the accuracy.


text:
gprofiler output ENSG to gene name
I would like to ask if there is any tool or R script that translates Ensembl IDs to gene names in a gprofiler GO output. An example of this file is the following one:

This is usually done using AnnotationDbi in combination with the annotation database for human org.Hs.eg.db.
library(org.Hs.eg.db)
genes <- c("ENSG00000102055", "ENSG00000232605", "ENSG00000136828",
           "ENSG00000159197", "ENSG00000201394", "ENSG00000278580",
           "ENSG00000213440", "ENSG00000166816", "ENSG00000275895",
           "ENSG00000105723")
select(org.Hs.eg.db, genes, c("ENSEMBL", "ALIAS"), "ENSEMBL")

This returns:
           ENSEMBL     ALIAS
1  ENSG00000102055       I-4
2  ENSG00000102055  PPP1R2P9
3  ENSG00000102055   PPP1R2C
4  ENSG00000232605      <NA>
5  ENSG00000136828   RALGEF2
6  ENSG00000136828  RALGPS1A
7  ENSG00000136828   RALGPS1
8  ENSG00000159197     ATFB4
9  ENSG00000159197      LQT5
10 ENSG00000159197      LQT6
11 ENSG00000159197     MIRP1
12 ENSG00000159197     KCNE2
13 ENSG00000201394   RN5S140
14 ENSG00000201394 RNA5SP140
15 ENSG00000278580      <NA>
16 ENSG00000213440      <NA>
17 ENSG00000166816     DLACD
18 ENSG00000166816       DLD
19 ENSG00000166816      LDHD
20 ENSG00000275895   U2AF1L5
21 ENSG00000105723     GSK3A

If you want to apply a 1:1 mapping of gene IDs to common names, you can use the "SYMBOL" column. In a data.frame, you can do something like this:
ensembl_to_symbol <- function(e) {
  hits <- select(org.Hs.eg.db, e, "SYMBOL", "ENSEMBL")$SYMBOL
  # Gene IDs without symbols should just stay as IDs
  no_hit <- is.na(hits)
  hits[no_hit] <- e[no_hit]
  hits
}

recode_symbols <- function(s) {
  e <- ensembl_to_symbol(unlist(strsplit(s, split = ",")))
  paste(e, collapse = ",")
}

df <- data.frame(
  genes = c(paste(genes[1:5], collapse = ","),
            paste(genes[6:10], collapse = ",")
  )
)

df$genes <- sapply(df$genes, recode_symbols)

This isn't the most vectorized (fast) solution, but it's probably easier to understand this way.

Yes I know the package. My problem is to use it in the format of gprofiler file. How could I  isolate each ensemble ID and replace it with a gene name, when ensembl ids are in one cell of a column of dataframe?THANKS!

I added an example of how you can use this in a data.frame.

Thnak you so much!!!


text:
Can you give me a detailed explanation about PrimeBlast?
I am trying to write a research paper on PRIMER-BLAST. I have read something about it on
https://www.ncbi.nlm.nih.gov/pubmed/22708584/
but still, it is not enough. So can you please tell me more about it or suggest some resources.

I'm not sure if it's just the language barrier and you actually meant something else, but...if you are unaware of basic tools for finding scientific papers (such as Google Scholar, Web of Science), and of good practices (such as following up citations of interest), then perhaps you shouldn't be writing a research paper just yet?

I would start at the Primer3 website here.
Once you understand all the variables on the website via the manual here and the associated publications Koressaar et al (2007) and Untergasser et al (2012) you will understand the gist of the algorithm used in primer design. This gives you optimised primers for a single strand of DNA/RNA
Primer-Blast looks in addition at maximising sequence conservation within an alignment of different sequences. The alignment is created using a Blast (I think its Mega-Blast. Thus Primer-Blast looks at one more parameter accounting for mutation heterogeneity which might be encountered in the PCR/RT-PCR and the number of primers it recommends is therefore reduced over say PRimer3.
Prove it for yourself take a given sequence and run it through Primer3 and then Primer-Blast and assess the differences and similarities.
References
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3 - new capabilities and interfaces. Nucleic Acids Research 40(15):e115
Koressaar T, Remm M (2007) Enhancements and modifications of primer design program Primer3 Bioinformatics 23(10):1289-91


text:
Model the effect of treatment on mutation burden
I am struggling a bit to model the following problem. Basically, I would like to model tumor mutation burden (mutations per megabase, a continuous) as a function of treatment (categorical). The following is an example of how data looks like
patient-id therapy-1 therapy-2 therapy-3 cancer-type TMB
PID-1      1         NA        1         CT-1        1.3
PID-2      0         0         1         CT-1        0.14
PID-3      NA        0         1         CT-2        5
PID-4      1.        1         NA        CT-3        2.3
PID-5      1         1         NA        CT-3        10
...
...

1 - treated; 0 - not treated; NA - No information
I was trying multiple linear regression, but in most cases one patient got multiple therapies. This is where I am confused on how to interpret the results of linear regression. So far, I haven't included the cancer-type information into the model. However, IMO, including cancer-type into model will yield specific treatment effects in specific cancer type? Or is it over-expectation of what we can achieve from linear model?
Also, is there a better way to address the problem to get reasonable associations between treatment and TMB in this case?
Thank you for the help. Stay safe!

Could you post your model is it something like TMB ~ cancer-type*patient-id ?

It seems you are trying to ascertain a relationship between treatment-type and tumor mutation burden (mutations per megabase). chi- square test for association could be helpful for a particular cancer-type. You can compute chi-square for each cancer-type.


text:
remove sequences from fasta file matching a string in the header
I have a file with 16S sequences. some headers contain species information. For my purposes I would like to exclude a number of species from the file, therefore I would like to do a pattern matching within a 16S header (that contains other characters beyond the species name) to then remove that header and the following sequence.
So, I have:
>NCBI_ID1;high_rank_taxonomy;...;species_name1
ATCCTT
>NCBI_ID2;high_rank_taxonomy;...;species_name2
ATTTTCCC
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA
>NCBI_ID4;high_rank_taxonomy;...;species_name4
CCCTTTAAA

I would like to exclude, say, species 1,2 and 4 to retain:
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA

I found this solution in an adjacent question:
awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

where I have my species names in the file ids.txt and sequences in file seq.fa
but that would work if the whole of my header was just the species name (contents of each line in ids.txt), which is not the case.
I thought I could be smart and use *species_name1* instead of species_name1 as contents of ids.txt but that apparently doesn't work.
Could someone help me with modifying the current solution to accommodate the existence of additinal strings in the headers?
Many thanks!
A

my awk is really bad.. would doing it in R or python help?

Say yes to that, @SW's Python is really good. The other way in is a Perl one-liner.

Or, you know, Heng Li's bioawk would be a good middle ground.

@StupidWolf sure thing, I am very new to bioinformatics, so I am happy to operate anything that does the job and learn it in the process : ) R would be preferable at the moment

Hi @Irrussional, see whether the R solution posted below works? I think it should be ok

yes, on it currently! thanks

you may use simple R commands, using the 'seqinr' package
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("filename.fa")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), ";"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c("species_name1", "species_name2", "species_name4")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fa")


text:
How to understand and analyse RNA-seq data (for a beginner)?
I am trying to understand expression of a certain protein across Pseudomonas species. I downloaded an SRA file from NCBI and converted it to a fastaq file. I am not able to understand how to interpret and analyse this file. I want to get the mRNA sequence and see which genes are being expressed and find the mRNA of interest.
How do I proceed with getting the mRNA sequences and then for mapping them to genome?

For some studies, in GEO they also include the processed counts. You can also see whether it's the case for your study


text:
What does instantaneous rate matrix mean?
PhyML is a tool used for maximum likelihood estimation analysis.
I was running a PhyML analysis to check the rate of nucleotide change in DNA sequences and in the PhyML output, there are rate parameters and an instantaneous rate matrix. Could you please clarify the difference between the two?

What I've understood so far is that the rate parameters indicate, for ex: that the A -C transition rate is 1.95 and so on. If so what does the instantaneous rate matrix tell us?

Hi @AntonyGuterres did you assign a root/outgroup to this tree before building it?

Reversible At its simplist the GTR rate, is the General Time Reversible model and most importantly infers the matrix is "reversible", thus as many A mutations will go to C as C mutations go to A and the number of mutations is relative to the G-T mutation rate. The G-T mutation rate is one of the least common mutations (transversion) so all other values are above 1. If you plotted this as a matrix the top right-hand matrix is identical to the bottom left hand matrix because its reversible, i.e. A <-> C . This is cool because you don't need to tell the algorithm which sequence is the root.
Directional The "instantaneous" matrix is a directional matrix, thus A -> C is not the same as C -> A. In order to obtain an accurate directional matrix you need to assign a root, its like lets start from the deepest node in the tree (the one between the root/outgroup and all ingroup sequences) and now - if that is correct - we can remove the restrictions of reversibility and impose directionality on the mutation process. This is needed if you are performing molecular clock type calculations. However, the key is assigning the root/outgroup - if it assigns it e.g. by mid-point rooting that is not necessarily a good idea.
Advice. My advice is if you are building a basic tree and haven't assigned a root/outgroup, obviously using GTR which is a good move, then you don't to worry about the intananeous matrix because it might not be meaningful. If you are running a molecular clock/dating/rate estimate type model and a very clear what the root/outgroup is and that is specified in the input, then the instananeous matrix becomes important.

The mutation rates of GTR and instananeous are not directly comparable because the GTR rate is with respect to G <-> T mutation, whilst instananeous is primarily concerned with the directionality of mutation and is a 12-parameter model, as opposed to GTR which 6-perameters, where one of those parameters is the reference.
Thus G-T is reversible in GTR (the reference) and so e.g. A<->C is reversible with respect to G<->T.
In instananeous matrix G -> T and T -> G are two different rates as are A-> C and C -> A, so one approximation based on an assumption with respect to a second approximation based on the same assumption can't directly equate to a mutation rate which removes that assumption.

I'm running a rate estimate model. I have supplied the guide tree for it to follow.
According to the Instantaneous matrix, A->C is 4.38 (and C->A is 0.153) but how is this related to the A<->C of 1.95 which is a GTR rate parameter?

They are not directly comparable answered above because A<->C (GTR) is with respect to G<->T mutation rate (which is assigned to a nominal value of 1).

Okay and how do you obtain A<->C (GTR) from A->C and C->A rates? Is it like a weighted average (weighted wrt base frequency)?

I think it is a separate question, I would do it via molecular clock trees.

Last two clarifications

In the matrix, N->C is high for all N. Does this mean that the a mutation most commonly results in 'C' vs any other base?
Also, a value at (x,y) in the instantaneous matrix represent the probability of base x mutating to base y right?

First question, yes, for the root you have set for the tree (it can change if you used a different root). Second question, yes and likely represent (row, column) - which is important in a directional matrix.


text:
Is there a software that can predict protein-DNA binding?
Is there any software that, given a DNA sequence as input, can predict what protein/ type of protein can bind to the DNA. The sequence of interest is a potential transcriptional regulator, I want to know for sure (or at least be reasonably certain) before I go to the wet lab to run any experiment.

try to do PROKKA predict protein.

Given a crystal structure or model of a protein it is incredibly hard predicting the DNA/RNA sequence it binds. If you do not know the protein, as is your case, homology based methods are the only way. And even then you may have a somewhat conserved non-coding region or motif, but nothing known because an analogue binds there. Transcription factors are present at very low levels so there is not a wealth of data to go by unfortunately, so your options are very limited.

@MatteoFerla Thanks for the clarification. What kind of homology based approaches can I try, are there any open source softwares

@AnthonyGuterres, what I meant to confer was that there is nothing as helpful as you had hoped. PROKKA will annotate a genome (if unannotated) and tell you what are transcription factors by homology to others. Some TF have known binding motifs (very organism dependent), so you'd filter for those with homology to one with a known motifs and go from there —unlikely to give you anything.
If it is bacteria, especially thermophilic, you can use guilt by association, i.e. TF is likely nearby on a different transcript. Works sometimes with funghi, while animals and plants are a mess.


text:
How to modify DNA evolution model to fit actual data?
I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).
I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms.
I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.
The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?

If mutation is a random process, why would you expect to reproduce the same result at the end of a stochastic simulation? Most simulation methods will not give you the exact same results as the training data; the "truth" is only one instantiation of many possibilities. Can you say anything about the methods you are using?

@MaximilianPress I agree, there will not be a 100% fit. Even 90% similarity is good enough.
Using a phylogenetic tree, every time there is a branch, I copy the genome and evolve it independently based on branch length. i.e. if a random number is less than branch length, the nucleotide changes (for each position)

What is the purpose of the sequence simulation?

@AnthonyGuterres see my new answer, it was going to be too hard to fit into a comment. Again, hard to know how to modify your model without knowing anything about your model.

Mutation is not a random process because there are many mutational biases, e.g. the transition vs. transversion ratio. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

INDELible http://abacus.gene.ucl.ac.uk/software/indelible/
Evolver https://www.drive5.com/evolver/ (not to be confused with the PAML package of the same name)

I decided that this was too involved for a comment, so I'm making it an answer.
It is still quite difficult to answer the question given the absence of methods explanation; but I think that there's some misunderstanding of phylogenetic modeling going on here, and I can at least address that. I'm going to try to break it down into its components, starting with the basic statistical idea of simulation and then going back to phylogenies.
Statistical reasoning
Let's say that I flip 10 coins (2 outcomes: head (H) and tails (T)), and I observe the sequence HTHTTTHHTT. I then fit a maximum likelihood model, and find that the best model for the data has
$Pr(H) = 0.4$,
$Pr(T) = 0.6$,
and that these two outcomes are mutually exhaustive:
$Pr(H) + Pr(T) = 1.0$
I can then use a binomial model to re-simulate data based on these parameters. I will observe various sequences, for example:
HTTTTTTHHT
TTTTHHTHHH
HTHHTHHTHH...
However, based purely on the parameters of the model, I don't expect every new simulation to recapitulate the exact same HTHTTTHHTT sequence, because I have reduced its information to some rate parameters for the simulation.
Each of those new sequences will yield model parameters ($Pr(H), Pr(T)$) that are similar to the training data, but they will not necessarily be particularly similar to the training data in sequence.
Phylogenetic modeling
Phylogenetic modeling works in the same way, such that you reduce the observations at the tips of the tree into some model with some parameters. Most of the time you are estimating e.g. a rate matrix of parameters for transition between states for a set of residues (in simplest case for DNA 4 nucleotide states and for proteins, 20 amino acid states usually). If you are working with gene sequences you are probably using the 4 nucleotide states or some derivation thereof.
Simulations usually work by starting with some fixed ancestral state at the root (which can take various states), and then evolving forward from there along each branch of the tree using the rate matrix in a continuous-time Markov chain to "replay the tape of life", by assuming that the model you've fit somewhat reflects truth.
While opinions differ somewhat on what to expect from "replaying the tape of life" (see that linked review for more), basically no one expects to get the same results as the input data, or even 90% similar to the input data. That would indicate that evolution is highly deterministic, which just isn't reliably true. Of course there are situations where such determinism happens to be the case, but they are relatively rare, and I don't think that sequence evolution is considered to be very deterministic at all.
It's a little different in that the tree defines a correlation structure of the tips, which constrains the simulation results to be somewhat more similar to the training data than we expect from coin-flipping. But there's still no reason to expect particularly similar results on the tips.
But you can also model discrete characters on a tree, which can be even simpler. For example, I used the R APE package to simulate the following character states (A or B for each tip) on a simulated tree of 10 tips with the following commands:
tr = rtree(10)
c_tr = rTraitDisc(tr)

I can then build a model for this character and this tree using the following:
rates = ace(c_tr, tr, type="discrete")

and use it to simulate new sets of states:
# A (index 1) is estimated to be ancestral state with 99% confidence
new_sim = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim, adj = c(1,0))

Look! Not a single observation of state B! Let's try it again:
new_sim2 = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim2, adj = c(1,0))

Ok, this is more similar, but there's still another state B on the tips.
Conclusion
Evolution is not very deterministic. The correlation structure of the tree makes the outcomes less random than independent trials as in coin-flipping, but it is not at all expected that simulation results from phylogenetic models will match the training observations well. The best that can be said is that each simulated dataset can hopefully be reduced to the same model parameters as the input data, but that's a very different statement from yielding the exact same outcomes.
If what you are looking for from simulation is to yield the exact same sequence states at the tips, then you are unlikely to get that from phylogenetic simulation. You might occasionally get pretty close, but that is not at all the expectation.


text:
How to annotate optimally a fungal genome without RNA-seq evidence?
Genome information：

~50M nt
2300+ contigs
No pre-trained parameters in Augustus
There are several well annotated RefSeq genomes of this genus.
There are several RNA-seq data of this genus, but not in this species.

My current strategy:
EST evidence:

I selected many RNA-seq data that include as many species as
possible, and de novo assembled these data using Trinity.
Only the longest transcripts were saved.
All transcript files were catenated into a Total_transcript.fasta file, and cd-hit-est was used to reduce data size with options -c 0.8 -n 5 -s 0.8.

Protein evidence

All fungal proteins from SwissProt was download and reduced using cd-hit with options -c 0.8 -n 5 -s 0.8.

Proteinsets from RefSeq genomes of this genus were download, catenated, reduced using same methods.

I planed to annotate genome using Maker2 with above Evidence.
I have a doubt about my strategy: Maker2 calls Augustus to ab-initio annotation, but no pre-trained model. How should I do to solve this?
Anyone else could share your valuable suggestions? Thanks anyway.

The Maker documentation does include some instructions for training ab initio gene predictors, but it assumes an abundant EST database is available. (Assumptions about what kinds of sequences will be available for a draft genome assembly have changed drastically in the last decade.)
It may be worth exploring an iterative approach in any case. You can do a first-pass annotation with Maker using the evidence you discussed, along with ab initio predictions from Augustus (and I'd recommend SNAP as well). Each gene model will be scored by Maker with the annotation edit distance (AED), representing how well the gene model agrees with the evidence (0.0 is best). From this first pass, if you can identify a few hundred reliable gene models with good AED scores (or maybe even as many as one or two thousand), that should be plenty of training data. You could use these gene models to train a new model for Augustus, using Maker's instructions or Augustus' own instruction. (In my experience, it was quite a bit of work in either case.)


text:
Mrbayes stop at certain generation without error
I am using Mrbayes v 3.2.7 MPI version to infer phylogeny, and I set the ngeneration to 10,000,000, but it often stop at certain generation (usually far smaller than given generations). Even if I restart the analysis by invoking the append=yes , but it stops again.
Anyone have this problem? I would appreciate you for sharing your useful experience.
Best regards,
Yanpeng
Version info

Version:   3.2.7
Features:  SSE AVX FMA Beagle readline
Host type: x86_64-conda-linux-gnu (CPU: x86_64)
Compiler:  gnu 7.5.0

Command line:
nohup mpirun -np 8 mb-mpi concatenate_phylipID.nex < /dev/null > /dev/null &

Mrbayes block in nexus file:
begin mrbayes;

    outgroup t1;

    charset Subset1 = 1-373 1953-2032;
    charset Subset2 = 374-531;
    charset Subset3 = 532-677 1583-1618;
    charset Subset4 = 678-1442;
    charset Subset5 = 1443-1582 1619-1725;
    charset Subset6 = 1726-1802 1803-1952;
    charset Subset7 = 2033-2661;

    partition PartitionFinder = 7:Subset1, Subset2, Subset3, Subset4, Subset5, Subset6, Subset7;
    set partition=PartitionFinder;

    lset applyto=(1) nst=6 rates=invgamma;
    lset applyto=(2) nst=6 rates=invgamma;
    prset applyto=(2) statefreqpr=fixed(equal);
    lset applyto=(3) nst=6 rates=invgamma;
    lset applyto=(4) nst=6 rates=invgamma;
    lset applyto=(5) nst=6 rates=invgamma;
    lset applyto=(6) nst=6 rates=gamma;
    lset applyto=(7) nst=6 rates=invgamma;

    prset applyto=(all) ratepr=variable;
    unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all) tratio=(all);
    mcmc nchains = 4 ngen = 10000000 samplefreq = 10000 printfreq = 10000 diagnfreq = 10000;
    sump;
    sumt;
    quit;
end;

How are you deciding it has stopped prior to 10,000,000 generations? Is it by counting the number of samples?

Check the log file,  It tells you how many generations have been done.

Have you tried running it in a screen or similar to see if it is something weird with nohup?

I prefer running prpgroms background if the analysis will take a long time. I think you are right. Ever I run it in PC five days, it works. Thanks.


text:
Number of sing-copy genes identified by orthofinder is far less than BUSCO assessment?
Dear Orthofinder users and BUSCO users,
Recently, I have assessed 30 genome annotataions using BUSCO with lineage dataset ascomycota_odb10, I got 1100+ complete and single-copy genes. But the orthofinder  only identified 300+ single-copy genes. Aneone else have encountered this similar problem?
Thanks,
Yanpeng

I know what caused the difference. BUSCO uses more strict criterion to determine single-copy genes but manly of the result genes are identified as multicopy genes by OrthoFinder using less strict parameters.

It's hard to say for sure without seeing the exact commands you used. However, Orthofinder identifies/defines single-copy orthologues as those genes present in a single copy in all species being compared in the analysis. If a BUSCO orthologue is present in more than one copy in one of the species being compared (could be multi-copy, could be fragmented etc.), then it won't be recovered as single copy by Orthofinder.

Thanks for your reply. I run OrthoFinder with default parameters. The BUSCO identified 1000+ complete single-copy genes but only 387 identified by OrthoFinder. I am trying to confirm whether BUSCO single-copy genes are true single-copy.


text:
Why all values become 1 after dcast?
I am reshaping my dataframe as follows.

I want to set the Wavelength column as the column names with the code as follows.
dcast(unique(df1, by=cols), ID + Condition ~ Wavelength, value.var = 'mean')

But all the values became 1s as follows.

I would appreciate any help.
Yanpeng

I think this should be posted on Stack Overflow - you’d get a better answer there anyway. But please make sure to post data as plain text and not as an image, otherwise it’ll likely get closed :) thanks

That comment is right, it's a general programming question, and text would be much easier for people to help with rather than images.
But, since I ran into this not long ago...
Take a look at the message shown in blue just above the output table.  My guess is you have some repeated ID/Condition/Wavelength combinations and it's defaulting to reporting the number of matching rows in the long-format data frame.  You'll need to either make sure you have all unique combinations, or explicitly define what you expect it to do with multiple values going into a single output row/column position (like, median/mean/max/min/whatever).  The unique() call won't do it since varying mean values would still make for unique input rows.
For example compare this:
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6222, 6223),
+     mean=c(106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
    ID Condition    662    664    666
1 6221      late     NA 106817     NA
2 6222      late     NA     NA 110357
3 6223      late 112998     NA     NA

with this:
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6221, 6222, 6223),
+     mean=c(106800, 106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
Aggregation function missing: defaulting to length
    ID Condition 662 664 666
1 6221      late   0   2   0
2 6222      late   0   0   1
3 6223      late   1   0   0


text:
no result from heat map WGCNA
I have created a heat map using WGCNA using R. But the heat map command is not showing the plot and there is no error. I have a large number of modules and traits. I am not sure if that is the reason. But I tried making it in HPC R and it still is not showing any result plot. Thank you!

Hi @Jessica you could post the code you have used, however the answers below do appear to have resolved the issue.

The plot is probably being drawn to a device that you cannot view on a cluster. If you're using the plotNetworkHeatmap function, say, try something like this:
png("plot.png", ...) #specify height, width, etc
plotNetworkHeatmap(....) #call your function
dev.off()

Now, scp the plot.png file over to your local machine to view the plot.

If you've copied over from some of the WGNCA tutorials you might have additionally copied the sizeGrWindow command - it's caused some of my plots to disappear.
ie:
sizeGrWindow(9,6)  # delete this line
plot(...)


text:
Is there a single cell RNAseq equivalent of GTEx or TCGA?
Or do I need to find individual studies and obtain data the long way!

Does the question read as "Is there a database of Human scRNA seq data?"? I think https://bgee.org/ contains some scRNA seq data, probably other databases too. Did you try to simply google it?

I think the most well-known scRNA dataset that is similar to GTEx in the Human Cell Atlas (https://www.humancellatlas.org). It is pretty much still ongoing but there are some preliminary data.


text:
How to filter a genome assembly consistsing of a large number of contigs?
I did some de novo genome assemblies with Illumina PE data using SPAdes, whereas most of them consisting of a large number of contigs(>1000). I have several questions below.

Do we need to filter short contigs?
How to determine the appropriate length threshold？

Hi, YP CHEN, could I know what kind of genome it is? bacterial? viral? eukaryotic? How long is the genome?

Fungal genomes， genome sizes are about ~ 50M

Thanks. I can't help in that case. But I guess what determines your length threshold is what use you intend to make from that assembly.


text:
How to query NCBI (Nucleotide database) by a feature qualifier?
I cannot seem to find a way to query nucleotide database (https://www.ncbi.nlm.nih.gov/nuccore) by a specific qualifier of the feature. For example, by /host under the feature "source".
For example, below is a part of one record, where the qualifier /host contains "Bemisia tabaci B".
How could I search for all sequences that contain the word "Bemisia" in the qualifier \host?
FEATURES             Location/Qualifiers
 source          1..351658
                 /organism="Candidatus Portiera aleyrodidarum BT-B-HRs"
                 /mol_type="genomic DNA"
                 /strain="BT-B"
                 /host="Bemisia tabaci B"
                 /db_xref="taxon:1206109"
 gene            1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
 CDS             1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_014894920.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="molecular chaperone DnaK"
                 /protein_id="WP_014894920.1"
                 /translation="MIRIIGIDLGTTNSCVAVLENGKSRVIENSEGARTTPSIIAYTD
                 DNEILIGQTAKRQAVTNPSNTLYAIKRLIGRKFNDDIVQKDIKMVPYKIVSAENGDAW
                 VQVKNRKLAPPQISAEILKKMKKTAEDYLGEKITEAVITVPAYFNDSQRQATKDAGRI
                 AGLNVKRIINEPTAAALAYGMDKNRGDKKIAVYDLGGGTFDISIIEIADVDGETQFEV
                 MATNGDTFLGGEDFDLKLINYIVKQFKLESGVDLSGDSLAMQRLKEAAEKGKIELSSS
                 QQTDINLPYITADKTGPKHLTLKITRAKLESLVEDLIKRSLLPCKKAIDDAKLTTKEI
                 DDIILVGGQTRMPLVQKKVAEFFGKEARKDVNPDEAVAMGAAIQGGVLGGDVKDVLLL
                 DVTPLTLGIETMGGIMTPLIEKNTTIPTKKTQIFSTAEDNQTAVTIHALQGERKKALQ
                 NKSLGRFDLTDIPPAPRGTPQIEVAFDLDANGILNITAKDKATGKEQSIIIKSSGGLS
                 EEEINKMIQDAETNKEEDKKFEELVQTRNQADGMIHTIRKTLKEIKNEDEKNKLESLI
                 NNLEQTLKSDDIKTIKEDLNKLTEEYSTIYQKIYSEKKPQNHNSEPAQKNAEKSSSKK
                 EGDVVDAEYEEVNNKTK"
                 ...

Best,
Milan

Doesn't it work if you just search for [host]"Bemisia tabaci B"?

@terdon Unfortunately it does not. If I try exactly as you wrote, it says "The following terms were ignored: [, ]".
If I reverse the query to be "Bemisia tabaci B"[host], then the response is "Unknown field was ignored: [host]."

Yes, it says The following terms were ignored: [, ], but the results seem to be what you're looking for, are they not?

Almost, it does include sequences that have the term in \host, but it also returns many sequences that do not have \host qualifier at all, but contain the term elsewhere, usually in the name of the sequence (DEFINITION). This is especially visible if you do [host]"Bemisia".

For this particular query, try host=Bemisia tabaci B[WORD].  However, you may also be interested in related queries like host=Bemisia tabaci B biotype[WORD].
But, I am not sure a general solution exists, as I do not believe all qualifiers are indexed.

Also, you can see how your query is interpreted on the right side under "Search details".  For the original attempt ([host]"Bemisia tabaci B"), it is interpreted as:
host[All Fields] AND "Bemisia tabaci B"[All Fields]

while the goal is to produce
host=Bemisia tabaci B biotype[WORD]

Thanks for confirming it to me, the issue is indeed that this qualifier is not indexed. In the end, I downloaded the results as GenBank file for a much wider search, and then locally imported into BioSQL database where I was able to write queries for every qualifier and get the result I needed. It takes some effort, but the flexibility of running complex queries was amazing to have.


text:
Gromacs molecular dynamic simulation throwing error during minimization
I am trying to investigate a membrane proteins influence on a lipid bilayer using molecular dynamics. I am new to MD and trying to run my first simulation. I have built the simulation input files using CHARMM-GUI and now I am using gromacs to perform the simulation.
CHARMM-GUI provides these commands to get gromacs to perform the initial minimization
gmx_d grompp -f step6.0_minimization.mdp -o step6.0_minimization.tpr -c step5_input.gro -r step5_input.gro -p topol.top
gmx_d mdrun -v -deffnm step6.0_minimization

The first command runs fine. The second command starts the minimization, but at step 15 it freezes.
...

Step 15, time 0.015 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms -nan, max 0.098209 (between atoms 16596 and 16597)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
  15683  15684   39.2    0.1109   0.1121      0.1111
  15683  15685   38.4    0.1115   0.1107      0.1111
  16596  16597   72.1    0.1111   0.1220      0.1111
  16640  16642   35.3    0.1116   0.1120      0.1111
  16643  16644   51.9    0.1108   0.1136      0.1111
  16643  16645   47.2    0.1114   0.1079      0.1111

In step 15, the rms is -nan, so I'm guessing that is the problem, and I am assuming that is coming from the difference between the "current" and "previous" is greater than constraint length like on this line here
16596  16597   72.1    0.1111   0.1220      0.1111

The problem is I don't have any intuition on how to fix it? The pdb file is a homology model from phyre without any modification other than adding a single chain id to all the residues.
I'm using a 1 to 1 POPC and DOPC in both membrane layers with a CHARMM36m force field.
The previous steps also throw LINCS warnings about bonds exceeding 30 degrees, however they produce a real value for rms.
I've counted on tutorials to get this far and would appreciate some suggestions or good resources I can use to continue my simulation.
Thanks!

Answer from @user1271772, converted from comment:
Try on this site: https://mattermodeling.stackexchange.com/questions/tagged/gromacs. There's an entire GROMACS tag there.
Answer from @vishal-kumar-sahu, converted from comment:
Try running/resuming using -cpi the simulation again, I got NaN error and rerun kept continued.


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

file = open("prats.fasta")

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

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

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

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

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

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

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


text:
Available phenotype data from GTEx
I am using gene expression data from the The Genotype-Tissue Expression (GTEx) project, available here: http://gtexportal.org. I would also like to use phenotype data for each sample (especially age, sex, race, body mass index, smoking history).
Access to such information is controlled by the NIH, so you need to apply for access.
Before I go through that probably lengthy process, I would like to know if it is worth it. Are any of the phenotype data age, sex, race, body mass index, smoking history available on a per-sample level, if you have access? I could not find this information anywhere.

If helpful for anyone else, I found that all and much more are available. See below 2 links.
https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study_id=phs000424.v8.p2&phv=169091&phd=3910&pha=&pht=2742&phvf=&phdf=&phaf=&phtf=&dssp=1&consent=&temp=1
https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno_variable_summaries/phs000424.v8.pht002742.v8.p2.GTEx_Subject_Phenotypes.var_report.xml


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


text:
How to download RNAseq gene expression data from GTEx
I would like to download RNAseq gene expression data from GTEx from the liver only. What is the best way to do this?
I have tried a few things. I downloaded several files from the GTEx website (https://www.gtexportal.org/home/datasets, under "RNA-Seq Data"), but the problems I run into are:

Not sure what is the best file: "Gene read counts" or "Gene TPMs"? What is the difference?
I would like 1 value per sample per gene. There are many duplicate genes in these 2 files. How can I collapse them to 1 value per gene?
How can I download the data for liver tissue only, instead of ALL data for ALL tissues?

Thank you for suggestions.

What do you mean by duplicates? I assume you mean multiple entries with the same gene name. Gene names are not unique identifiers. Also, whether you want raw counts or TPMs depends on what you want to do with the data. If you want to calculate differential expression between 2 samples, you will need counts. Honestly, comparing expression values of different genes is apples to oranges. Plus, you need replicates, ideally from the same batch. RNA-Seq is very noisy data. Finally, if you want a specific tissue, you have to find samples of that tissue.

@ATpoint and hepcat72  Thank you. I have found that subsetting the data worked to get data from one tissue only. The gene names start with ENSGxxxx. I found out that they are different transcripts, sometimes of the same gene, so I suppose there are no real duplicates.

Gene read counts are the number of reads or fragments that are mapped to each gene. So it is an integer and the values are not directly comparable between samples due to differences in sequencing depth. TPM stands transcripts per million. It is the normalized gene expression level. Basically, you first normalize read counts by gene length and then normalize the library size to 1e6 in all libraries. As for which one should you use, it depends on what analysis you want to do and what tool you are going to use.
I have not noticed that there are duplicated genes in those files, can you show an example?
I do not think you can download expression from one tissue, but you can always subset. The whole dataset is not very large.

Thank you very much. Subsetting worked for me.


text:
How to get cytoband and gene level copy number from genome wide SNP array copy number data?
I have (human) Illumina genome wide SNP array copy number data. For each SNP genomewide, I have Log R Ratio (LRR) and B Allele Frequency (BAF).
What tool(s) can I use to get the integer copy numbers (either -3 to +3 or 0 to inf) for each cytoband and for each gene in the genome?
I assume that people routinely do this? Are there any R packages/pipelines/code to do this?
Do people use DNAcopy and/or GISTIC for this?
Thank you.

To obtain cytoband and gene-level copy numbers from genome-wide SNP array copy number data, you might want to consider using the ASCAT tool, as it makes allele-specific copy number segment data and derived integer gene-level copy number values from SNP array data. The SNP6 copy number analysis pipeline (ASCAT2) I believe is specifically for this purpose.
See also:

CNVkit, which has a method for copy number detection via using targeted reads and the non-specifically captured off-target reads to infer copy number evenly across the genome.

Array Comparative Genomic Hybridization (aCGH) and SNP Genotyping Technologies

CNVannotator

Thank you very much for your answer with multiple options. I will look into these. ¡Gracias! :)


text:
Interpretation of clusterProfiler results
I am using clusterProfiler for differential gene expression using this code:
https://yulab-smu.top/biomedical-knowledge-mining-book/enrichplot.html
I have a question about the third line of code:
 de <- names(geneList)[abs(geneList) > 2]

Here all gene names are taken with absolute fold change >2. As an example, if in a hypothetical extreme case you have 100 downregulated genes (log2 fold change <-2) and zero upregulated genes (log2 fold change >2), then you may incorrectly think that identified pathways are activated, but they are actually repressed. How is this accounted for? Thank you.

Hypothetical scenarios don't make the best questions, but this question's ok. Do keep in mind this is differential gene expression rather than explicitly being repression (don't always intersect).

I agree with you :)

The scenario you described would be reported by the code because:
... [abs(geneList) > 2]

The abs will change any negative value to positive and retain all positive values as positive. Therefore values <-2 will be reported as well as values >2. abs = absolute and to take the absolute value of a number is exactly the operation just described.

From comments. Just remove abs
negative <- names(geneList)[geneList < -2]
positive <- names(geneList)[geneList > 2]

Thank you for your feedback M__. Yes, this operation selects all genes with >2 and all genes with <-2 fold change in expression. However, none of the values are themselves retained. It is only the names of the genes that are taken forward to the next step. Therefore, the direction (- or +) is lost after this step. In other words, if you have gene A +8 fold change, gene B +3 and gene C +2.4, then the end result is the same if gene A had been -6.2 fold change, gene B -2.4 and gene C +2.1. I suppose I conclude that we are looking at representation, rather than at up/downregulated pathways?

@SylviaRodriguez addressed in the comments

Thank you. I like that solution. Two separate analyses.

Thanks @SylviaRodriguez I'm glad its resolved.


text:
What is the difference between fixed effects and random effects in the context of Linear-mixed models?
The terms in LME i.e. fixed and random create confusion? What is the genesis that can distinguish between the two ?

A dog's response to medicine may depend  on number of days it is sick and  its sex or breed.

LME assumes two components- fixed and random. Fixed effect (mean) is constant and consistent. THE random part shows a number of variables that influence the dependent effect-size individually.

In case, we have predictor variables of continuous as well as nominal(or level-type 2 or 3) categories data, it becomes necessary to use linear as well  as non-linear (probabilistic) models. Therefore, Mixed models convey use of two types of models - linear and non-linear. Nonlinear models work under probability and linear models perform under fixed effect plus random effects.  The goof up is Mixed model.

https://stats.stackexchange.com/questions/4700/what-is-the-difference-between-fixed-effect-random-effect-and-mixed-effect-mode

Maybe worth reading the link and the article above.. I have no idea where this question is going without a context, and what is its relevance to bioinformatics?

This is a good accessible description and explanation of what random effects are https://gkhajduk.github.io/2017-03-09-mixed-models/

The definitions e.g. fixed effect  and random effects are not congruent with statistics and mathematics. Any way, thanks for posting the link.

https://www.stat.ubc.ca › BookPDF
Web results
Mixed Effects Models for Complex Data - UBC Department of Statistics - University of British ...

I do not know much about statistics but I will try my best to explain.
First, random effects are defined as the factors (categories) in the population that we are not aware of (not observed), so we are randomly sampling levels of those factors when we sample the population.
Practically speaking, random effects can be found when there are hierarchical structures in the data. For example, we want to compare the relationship between students' SAT score and their college admission rates. You might just want to run a simple linear regression. But there might be some structure in the data that is hidden. For example, the students can be grouped into classes, and then schools. You can also group students by the social-economic status of their parents.
The idea is that students belonging to different categories may have different college admission rates even if they have the same SAT score. Therefore, it is important to consider some of these "random effects" when building your model.

Linear regression hypothesis is that  effect of a
 variable is comprised as a constant component plus variability. In other  words,  the independent variable (s) is deemed to reflect a fixed (constant) part - determined as population effect size - and a variable part - random effect. Further, the term -
Linear Mixed(regression)model  is used to reflect use of continuous variables and nominal variables acting as independent variables that influence  the dependent  variable.  THE CONTINUOUS VARIABLE  is  comprised  of population effect size - and a variable part - random effect . The nominal or categories based variable effect  is probabilistic and we use the term probabilistic variable for such variables. Sometimes academics mistakenly call its effect on dependent  variable as random effect.
Finally, this may be noted that  principles of ANOVA and probability  have been combined to formulate linear Mixed effects models.


text:
Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper
I am trying to replicate the results of this paper roughly guided by this pipeline.
Basically, we are trying to detect a differential expression between mesenchymal and epithelial cells under cisplatin treatment vs control.
However, I am facing a problem when trying to calculate the p-value for the Mann-Whitney U-test. This is the part of the paper I am referring to:

Robust Multichip Average (RMA) normalization was performed at the
transcript level on the results from the Affymetrix Human Gene ST1.0
arrays using Affymetrix Power Tool 1.14.4.1 for all 46 sham- or
cisplatin-treated ovarian cancer cell lines. The normalized data were
subsequently standardized using ComBat71 to remove the batch effect.
In this experiment, the cisplatin treatment assay was performed in
triplicate on 20 cell lines, while single assays (without replicate)
were performed on the remaining 26 cell lines. Taking advantage of the
triplicate data, potentially fragile probes with strong variations (an
s.d. of >0.2) within the triplicates were removed, decreasing the
probe number from 33 297 to 21 329. To perform a fair comparison, the
triplicate data were then log-averaged into one value so that one
result for each cell line could be used in the following analyses. To
detect differential responses to cisplatin between epithelial- and
mesenchymal-like cell lines, the transcriptomic responses to cisplatin
were computed by subtracting the gene expression value of control
(cisplatin untreated) cells from that of cisplatin-treated cells.
Mann–Whitney U-test (P<0.01 as a cut-off value) was subsequently used
to detect the differential transcriptomic responses between the
expression changes by cisplatin treatment in epithelial-like cell
lines with those in mesenchymal-like cell lines (Supplementary Table
3).

My p-values for a given gene are astronomically different to that of the paper...
Here is my code:
expr_f is the expression values matrix from the Expression set (log-transformed).
SDRF_f is the pData from the same expression set.
I am feeding into the wilcox.test function a vector of expression from epithelial cells (epi_values) and mesencymal cells (mes_values) for the same probe (8143663).
 genelist=rownames(expr_f)
epi_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Epithelial-like"))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Mesenchymal-like"))

epi_values = expr_f["8143663", epi_samples]
mes_values = expr_f["8143663", mes_samples]

wilcox.test(mes_values, epi_values, exact=TRUE)

In the paper, their p-value for this probe is 0.000372259258165627 while mine is 0.8785. The difference makes no sense at all to me.
I don't know if I am not quite getting the concept of the u-test, if I am using the function wrong, or if it's something else. Any input is appreciated

Cross-posted on biostars: https://www.biostars.org/p/467695/

Wow that is some difference particularly for the Mann-Whitney. You are using Wilcox of course. Regardless there is one possibility which I will mention below if I get chance.

I think you may need to check your pipeline in this regard. I ran the analysis and replicated as closely as I could their M&M (I didn't filter the "fragile probes" since you anyway care only about a specific probe). In any case, I get something relatively close to their result (P=0.004059), which is still one order of magnitude off, but without more detailed M&M, this is probably as close as I can get to their analysis. Here's what I did (if you prefer to read the notebook html, that is here)
library(Biobase)
library(GEOquery)
library(sva) # ComBat
library(oligo)
library(googlesheets4)
library(stringr)
library(dplyr)

#' Get the GEO dataset mostly for the phenoData annotation
gset <- getGEO("GSE47856", GSEMatrix =TRUE, getGPL=FALSE)
adf <- Biobase::phenoData(gset[[1]])
sampleNames(adf)

#' Read RAW cel files from the GEO data supplementary rar archive
#' https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE47856&format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles <- sapply(sampleNames(adf), function(f) {
  pat <- paste0(f, ".*.CEL.gz$")
  dir("~/Downloads/cel", pattern = pat, full.names = TRUE)
})
raw_data <- oligo::read.celfiles(celfiles,
                                 phenoData = adf,
                                 experimentData = experimentData(gset[[1]]))

#' Get the raw expression and RMA on transcript clusters
rma_expr <- oligo::rma(raw_data, target = "core")
rma_expr_assay <- assayData(rma_expr)$exprs

#' Use ComBat for batch effect removal
batch_data <- phenoData(rma_expr)[["batch:ch1"]]
cb_expr <- sva::ComBat(rma_expr_assay, batch_data)

#' Get supplementary table 1 from the paper (which I copied on this Google Sheet)
supp_data <- read_sheet("https://docs.google.com/spreadsheets/d/142__1O35zvlttejvz62XQX33uLRJhqIDvb604smVSRo")

#' Match the 'cell line' phenotype from the phenoData to the sample info of
#' the supplementary table and extract the cell classification. Also remove
#' all white space and convert everything to lower case for easier matching.
#' Some cell lines also have a/b attached, so I'm dropping that to match them
#' up. Mayb not be correct to do that here and treat all a/b the same (?)
cell_line <- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[["cell line:ch1"]]),
    " +"),
  "^tay$", "taya"),
  "dov13[ab]", "dov13"
)

sample_info <- str_remove(str_to_lower(supp_data$Sample), " +")
classification <- supp_data$Classification[match(cell_line, sample_info)]

#' Get the treament from the phenoData object
treatment <- adf[["treatment:ch1"]]

#' We'll need the geometric mean later (R has no built-in)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

#' Set up a tibble (==data frame) with all that info for your probe
tibble(Expression = cb_expr["8143663", ],
       Line = cell_line,
       Treatment = treatment,
       Classification = classification) %>%
  # Summarise replicated data with the geometric mean
  group_by(Line, Treatment, Classification) %>%
  summarise(Emean = gm_mean(Expression)) %>%
  # Make sure data is in order
  arrange(Line, Treatment, Classification)%>%
  # Calculate the difference between the cell types
  group_by(Line, Classification) %>%
  summarise(EmeanDiff = diff(Emean)) -> xdata

#' Test the difference
xdata <- split(xdata$EmeanDiff, xdata$Classification)
wilcox.test(x = xdata$`Epithelial-like`, y = xdata$`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
```

Hi @CSP it is certainly an interesting discussion and some good R coding. Could you look through this very detailed reply and comment please?

You were right, I had a minor mistake in my pipeline! Thank you for this answer, very detailed and useful!

Great that it worked. If this is answered, please accept the answer

Agreed @CSP please accept the above answer

The most likely explanation here is paired versus unpaired.
Mann-Whitney is normally associated with unpaired data, but the authors may have performed a paired test. You performed an unpaired test and thus there is a massive discrepancy in results. The authors should have stated the effect of pairing - it that is what they did - but they are under no strict obligations to do so. You however have performed a Wilcoxan rank sign test, albeit the results shouldn't be so different.

In R, a two-sample wilcox.test (i.e. providing both x and y) is the same as a Mann-Whitney U test. From the docs:

Performs one- and two-sample Wilcoxon tests on vectors of data; the latter is also known as ‘Mann-Whitney’ test.

Also they cannot do a paired test, since they have uneven groups (Epithelial vs Mesenchymal)

Okay fair point @BastianSchiffthaler , I'm wrong. I was aware of the ability of Wilcox to approximate Mann-Whitney ... if the samples are naturally unpaired ... again simply I'm wrong. Mann-Whitney is naturally unpaired so I was never really sure. The only thing I'd point out is that generally any bug can result in p<<0.01 but when p>>0.1 it difficult for a bug to achieve this.

The samples are unpaired, it was a mistake in my prior pipeline, but thank you for you very useful input!


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

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

A FASTA file with only one sequence.

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

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

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

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

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


text:
Looking for a tool to find 16S rRNA in hundreds of genomes
I have hundreds of bacterial genomes and I would like to know which is the best informatic tool in which I can run all genomes on one go and get a 16S rRNA gene sequence of each genome as output.
Could anyone give me some advice?

Have you tried anything any particular? As a general rule, there is no single "best" tool for much of anything, you just need something to give you a reliable result.

You can use Barrnap
barrnap -k bac --threads "$(nproc)" \
        -o output_rrna.fna < input_file.fna \
        > output_rrna.gff 2> /dev/null

Put the file in one line of sequence per fasta id and get sequence with 16S annotation.
perl -pe 'if(/\>/){s/\n/\t/}; s/\n//; s/\>/\n\>/' output_rrna.fa | \
    perl -pe 's/\t/\n/' | tail -n+2 | cut -d\: -f1 | \
    grep -A1 16S > output_16S.fasta

Check the output_rrna.gff and make a BLAST of the sequence to be sure about the predicted 16S.

Thanks zorbax, I am going to have a look and I will give you a comment about what I get

You can perform genome annotation using tools like prokka or PGAP. Using some python scripts  you can extract the predicted 16S and blast to confirm.


text:
Tools for modelling and visualizing growth of cells
I am looking for a program/GUI app/package that would help me do simulations for cell growth. I have a microscopic mathematical model, and the scenario is basically the following: I start with a rod-shaped cell characterized my its center and direction angle, it grows until it reaches a certain length, it divides into two identical cells whose angles change randomly and the process continues.
Is there any program or package that I can use for such a situation? I have found several on GitHub, but they were mostly for very complicated biological processes. I use Python/Julia/Matlab and can code, but I was wondering whether there are any tools for such simpler scenarios.

If the scenario is simpler, and you can't find anything by searching, I would consider just writing your own.
The modelling can be done in any language but if compute time is an issue - for example you want to run a simulation with many cells or many simulations - then I would use a language that can support that such as Julia or NumPy in Python. In fact Julia is a nice language for this type of thing as you end up with readable code that runs fast.
For visualisation in Julia specifically you could use the animation functionality of Plots.jl, which is basic but there, or could try Makie.jl. See the Molly.jl documentation for some examples of visualising particle-type processes using Plots.jl.


text:
Visualizing cell growth
I model the following events: a rod-shaped cell with center $x$ grows symmetrically until it reaches a maximum length/age, it divides into two identical cells, and the process continues. At division, the direction angle $\varphi$ of the daughter cells changes randomly.
So, the growth rate is $l=\frac{l_{max}}{2}+\frac{a}{a_{max}}\frac{l_{max}}{2}$
Is there any way of naively visualizing this? I found different cell modellers on GitHub but they were very complicated biologically speaking. I have some basic Python knowledge, but I have never done simulations.

Are you sure it is not 1/l ?

I would recommend looking into Python's plotting utilities, for example here.
I would recommend specifically looking at the examples that plot mathematical functions.
The basic idea is that you create a function that computes your equation based on the different parameters as inputs ($a, a_{max}, l_{max}$) and then you can plug in a range of values for each parameter and plot e.g. with $a_{max}$ on the x axis, then different lines showing a range of values for $l_{max}$, or whatever makes sense for you.
e.g.
import numpy as np
import matplotlib.pyplot as plt

def fn(a, amax, lmax):
  return lmax/2 + (a/amax) * (lmax/2)

plt.figure()
amax = np.arange(0.01, 1, .01)
plt.plot(amax, fn(.005, amax, .1), "ro") ## red points
plt.plot(amax, fn(.005, amax, .2), "go") ## green point
plt.show()

I have no idea if these parameter ranges are meaningful for you, you will have to figure that out for yourself.

Thank you for your answer! I, however, messed up in the first place. The growth rate is obviously a velocity, equal to $$s=\frac{l_{max}}{2a_{max}}$$.
After doing some reading, the length of division should be fixed, and the cells are theoretically supposed to double every 30 minutes. However, as it gets crowded, the maximum age might become bigger.

Ok, I'm not sure that I understand what you are saying but you should be able to use the tools in the same way regardless of whatever the function turns out to be.


text:
NCBI SRA database sample: control vs test
I was trying to download some data from NCBI SRA (SRA059451). There are 27 samples available for SRA059451. But i am unable to understand which samples are 'control' and 'test' samples. Please help me regarding how to differentiate the samples.

Seems to be cross-posted on Biostars, which is fine but please refer both posts to each other so that answers on one site would be utilized by the other community.

You can download the information using this link:
I did it in R:
tab = read.delim("https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=SRA059451&result=read_run&fields=run_accession,experiment_title,sra_ftp,sample_alias")

head(tab)
  run_accession
1     SRR585697
2     SRR585700
3     SRR585726
4     SRR585727
5     SRR585728
6     SRR585729
                                                                experiment_title
1 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 171
2 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 238
3 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 241
4 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 242
5 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 243
6 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 245
                                      sra_ftp
1 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585697
2 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585700
3 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585726
4 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585727
5 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585728
6 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585729
                                    sample_alias
1     Symptomatic atherosclerosis gut metagenome
2 Symptomatic atherosclerosis gut metagenome 238
3 Symptomatic atherosclerosis gut metagenome 241
4 Symptomatic atherosclerosis gut metagenome 242
5 Symptomatic atherosclerosis gut metagenome 243
6 Symptomatic atherosclerosis gut metagenome 245

Normally the information will be here, but I cannot find what does metagenome 238 etc refers to. No table in the paper that can link sample to phenotype..
Most likely you need to contact the authors. This list might be useful though.


text:
SRA run data access
Trying to download run data of Run accession SRR2155174. Here's the link: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR2155174
Please go to the link and go to the "Data Access" tab. Hereunder header "Name" you can find two links:
https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR2155174/SRR2155174.1
https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-6/SRR2155174/SRR2155174.1
Do these two links download the same data or they are for different purposes?
Thanks

It seems just like two different mirrors contains the same SRA archive. If you still doubt, you can download both (its have relatively small size) and compare their md5sum.


text:
How to generate Beta diversity boxplot from phyloseq object?
I am working with profiled metagenomic taxonomy abundance data. I want to generate beta diversity (bray-curtis) boxplot from a phyloseq object where two groups (control and test) will be shown. Something like this figure:

How can I make it? Can anyone share any code or tutorial?
Thanks,
dpc

You could convert the phyloseq object to a dataframe and just plot it using base R or ggplot2. I'm not exactly sure what the bray curtis df will be called but you can access it using the @ notation then convert with data.frame e.g.
bray_for_plot <- data.frame(phylo_obj@bray)


text:
How to decide between standardized mean difference (SMD) and log2fold change (log2FC) before doing meta-analysis?
I want to do meta-analysis of alpha diversity between test and control groups across 10 different populations. I want to see what is the summed effect (pool effect) of intervention between the two groups for alpha diversity. For that, first I need the intervention effect on test group's alpha diversity for individual population dataset. How should I calculate this? Should I calculate this by difference in the log transformed mean values of the groups (log2FC) or should I calculate SMD for each study?
Thanks

Cross-posted on biostars: https://www.biostars.org/p/9468060/

what is the measure of alpha-diversity? counts or proportion?


text:
When should I apply a multiple comparison adjustment?
I am performing a RQL analysis ("http://www.esapubs.org/archive/ecol/E095/002/suppl-1.pdf").
This analysis meassures the relationship between different traits of three different matrix showing if there is a positive or negative relationship. This RQL analysis make a fourth-corner method.
Well, this method have the option to  make a p-adjust method to solve multiple comparison problem. However, my variables are independent and I belive that I should not apply to my dataset. Am I right? Or always that two or more matrix are compared I should apply a multiple comparison adjustment?
Thanks in advance.

The p-value adjustments corrects the inflated number of significant results when performing many significance tests.
If you have a test with 5 % false positive rate and but you run it several times the chances of getting a false positive results in any of these tests is far greater than 5 %. For example repeating the test with 5 % false positive rate 20 times means that there is a 60 % chance of at least one false positive result.
You can see this in the tutorial you have linked to in your question, where the number of significant differences goes down from 26 without correction to 18 with correction.
Therefore, if you are performing multiple comparisons you should do a correction. The problem of inflated p-values will occur when you are performing many tests even if the tested variables are independent.


text:
Error after trimming illumina adapters
I am removing illumina adapters of the NGS data with a loop. My NGS data is storage in /data/HTS_seq/.
I used this function:
 for infile in /data/HTS_seq/*_read1.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

 for infile in /data/HTS_seq/*_read2.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

and It gives to me this error:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'
  Traceback (most recent call last):
  File "/usr/local/bin/cutadapt", line 10, in <module>
cutadapt.main()
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 652, in main
options.discard_untrimmed)
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 327, in trimmed_and_untrimmed_files
trimmed = untrimmed = xopen(output_path, 'w')
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/xopen.py", line 178, in xopen
return open(filename, mode)
      IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which is my error?
Thank you in advance.

First of all, you aren't actually setting the $outfile variable. Your loop has this basic format:
for var in *; do newVar="something" command $newVar; done

For example:
$ for var in 1 2 3 4; do newVar="new:$var" echo "$var: $newVar"; done
1:
2:
3:
4:

As you can see, newVar is never actually being set. That's because setting a variable is its own operation. So you need to do that, end it with a control character like ;, or &, or a newline, and then run your cutadapt:
$ for var in 1 2 3 4; do newVar="new:$var"; echo "$var: $newVar"; done
1: new:1
2: new:2
3: new:3
4: new:4

Or, with your data:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

But that would give a weird file name, like foo_S1_L001_read1.fastq_trim.fastq. You probably want names like foo_S1_L001_read1_trim.fastq instead, so you would want:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Now, the specific error you show makes me think that isn't the command you ran. Could you perhaps have run this one instead by mistake (note the / instead of \ before the _trim.fastq)?
for infile in /data/HTS_seq/*_read1.fastq; do
 outfile=$infile/_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Your error message is:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which means that cutadapt tried to open /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/ as a directory and then look for a file named _trim.fastq inside it. The only way I can think that this would happen is if you ran the command I gave above and not the one you show in your question.

It works! Thank you Terdon

There is a bug here
/data/HTS_seq/*_read1.fastq

This will go to your root directory / and there is definiately no folder called data in the root. To prove it type ls / and that is where the command is pointing, then compare this to ls ./.
To solve it go to the directory HTS_seq and type pwd then copy and paste the directory structure. I suspect ./data/HTS_seq will solve it
 for infile in ./data/HTS_seq; ...

This will look for the folder data in your current working directory rather than the root directory.

Okay you've been assigned /data as root ... no biggie.
PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq

Why is / in the file name? What happens if you ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq?
and ls /data/HTS_seq/ Finally could you please ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq\/_trim.fastq?
The \/ is important. What I think has happened is a / is in the file name and that will cause problems because it will think that _trim.fastq is a directory when its part of a file name, and would explain the error.
Then do the same thing using head instead of ls see what happens.
If the last ls and head works then rename the files by typing rename in the command line. If rename works it will be easy to remove / from all your files using the power of the rename function (its a Perl function BTW).

Hi Michael, thanks you for your answer. I typed pwd in HTS_seq folder and I got this direction: /data/HTS_seq                                                                                Anyways I typed the code with the . :                                                                                for infile in ./data/HTS_seq/*_read1.fastq; do outfile=$infile_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done and I got the same error. Which is my mistake? Thanks you

Hi @AdriánP.L. please try the above, I should have read the error more carefully first time around

You can never have a / in a file name. The / alone with \0 are the only two characters that can never be in file names.