text:
calculateQCMetrics defunct, how to calculate the quality metrics by perCellQCMetrics
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use calculateQCMetric function to calculate the quality metrics, but I am getting an error message that calculateQCMetric is a deprecated (isSpike function is also deprecated). I want to calculate the quality metrics. Can anyone provide code for getting QCMetric step? I attached the previous defunct code below:
umi <- calculateQCMetrics(
   umi,
   feature_controls = list(
       ERCC = isSpike(umi, "ERCC"),
        MT = isSpike(umi, "MT")
   )

And I use StupidWolf's code from this question isSpike function in SingleCellExperiment package is deprecated? to instead isSpike function. I also attached below:
Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]
Have a vector of genes to keep:

keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))
Then make the single cell object:

sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)
Use altExp to slot in the stuff:

altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

I try the below code, which did not work.
> perCellQCMetrics(sce,
 subsets = list(Spikein_names,MTgene_names),
 flatten = TRUE,
 exprs_values = "counts"，
 use_altexps = TRUE)

Thanks in advance.

Isn't that basically the workflow from Scater? I would dive into the vignette and see how they perform the lowlevel QC.
https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html

I'm also following the Sanger tutorial. I've tried the code suggested by mattcwh but get the following error: Error in assay(altExp(x, use_altexps[i]), exprs_values) : 'assay(<SummarizedExperiment>, i="character", ...)' invalid subscript 'i' 'counts' not in names(assays(<SummarizedExperiment>)) It might have something to do with logcounts, at least based on this post (?) : https://github.com/LTLA/scater/issues/3

Ok I have figured it out and I am sharing this in case anyone can get help. I used
sce <-addPerCellQC(sce)

which calculated the QC for all the altexps also and adds them to the colData.

It would be helpful if you can provide the error message.
Guessing from your codes, it's possible that the alternative experiments weren't set correctly. My suggestion:
altExp(sce, "spike-in") <- SummarizedExperiment(assays=list(counts=as.matrix(SpikeIn)))
altExp(sce, "MTgenes") <- SummarizedExperiment(assays=list(counts=as.matrix(MTgenes)))

If you're using altExp, the subsets argument in perCellQCMetrics is unnecessary in this case.
perCellQCMetrics(sce,
exprs_values="counts",
use_altexps=TRUE)


text:
Why use Needleman-Wunsch if there is no way to evaluate the statistical significance
I thought that Needleman-Wunsch is the best approach to align sequences. However, I read that it is impossible to evaluate the statistical significance of the alignment if you do global alignment. So I would like to ask, what is the point of global alignment if we cannot evaluate statistical significance?

Can you provide more context to this? Statistical significance doesn't normally have anything directly to do with an alignment. Dynamic programming has nothing to do with any of this other than it's a convenient way of implementing things like Needleman-Wunsch.

I am sorry. In my notes, it is mentioned that in the case of global alignment you cannot know the statistical significance of the alignment, whereas in the local alignment without gaps you know that the distribution of score is Gumbel and you can calculate p-value, bit score and e-value. And yes when I said dynamic programming, I meant Needleman-Wunsch algorithm.

Can you say more about how these two things are related? We use dynamic programming in alignment just to speed things up I think. I am not sure whether people are interested in a p-value for global or local alignment.

i mean sequence alignment like in blast. to be sure that the alignment is not by chance

It's best to flip the question on its head and ask yourself instead, "Do we care about the statistical significance of an individual alignment?" The answer is "yes" if you're in a context like homology modeling, where you're aligning two protein sequences. The answer, however, is a very definite "no" in the context of NGS, where you just want the most likely alignment and a vague assessment of its likelihood of being wrong (for filtering purposes).

thank you!but i am interested in homology modeling and not in ngs. So alignment in order to find some relationship between sequences.

source: Stack Exchange

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

sort: votes


text:
Edit FASTA header using sed
I need to rename the following headers from a FASTA file. Something like:
>ENST00000488147.1::chr1:29570-30013(-)
>ENST00000405392.5::chr2:24919186-24920186(-)

into different format like this:
>chr1:29570-30013(-)_ENST00000488147.1
>chr2:24919186-24920186(-)_ENST00000405392.5

How can I rename the FASTA headers, perhaps using sed or awk?

Is there fasta sequence between the headers? If so is it multiline or single line fasta? Also what about Biopython?

Hi Michael. Yes, there are fasta sequences between the headers, and it has a large number of fasta sequences alongside their headers for each. I didn't know anything about Biopython, but I just looked it up and it seems really helpful as it has motif analysis tools (MEME or MAST). Thank you

awk -F "::" '{if($1~">"){gsub(">","");print ">"$2"_"$1} else {print $0}}' foo.fa

Basically the same as from Devon but using -F to indicate the initial field delimiter and then using gsub to remove the >.

Hello ATpoint. Thank you for your answer. I should've used gsub instead of sed to remove the string that I don't want. Thank you for your answer

awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa

An explanation of the non-obvious bits:

BEGIN{FS="::"} split columns using :: as the designator
$1~">" If there's a > in the line
substr($1,2) Trim the > that was at the beginning of the line.

awk -F:: '{sub(/>/, ">" $2 "_"); print}' KDM4A_1kbp_TSS.fasta > KDM4A_1kbp_TSS_swp.fasta         sed 's/::chr.*)//' KDM4A_1kbp_TSS_swp.fasta > KDM4A_1kbp_TSS_MEMEinput.fasta I just checked your answer now. I had tried working by myself and come up with this script. But, yours looks much better than mine :) I appreciate it for your answer.

sed -E 's/>(.+)::(.+)/>\2_\1/' in.fa > out.fa
If you want to use sed.

This misses the > at the beginning of each name:
sed -r 's/>(.+)::(.+)/>\2_\1/' for linux and
sed -E 's/>(.+)::(.+)/>\2_\1/'  for macOS.

Sorry there was a little typo. Now it should work properly

I would use perl here instead of sed so you can use non-greedy patterns (e.g. .*?) and so ensure that you always match the first occurrence of :: if there are more than one on the line. Perl also has -i, and in fact is where sed got the idea from, so you can edit the file in place just like you can with sed.
Using this example file:
$ cat file.fa
>ENST00000488147.1::chr1:29570-30013(-)
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>ENST00000405392.5::chr2:24919186-24920186(-)
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

You can do:
$ perl -pe 's/>(.+?)::(.*)/>$2_$1/' file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

Or, to edit the file in place:
$ perl -i -pe 's/>(.+?)::(.*)/>$2_$1/' file.fa
$ cat file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

Explanation

perl -pe: print each line of the input file after applying the script given by -e to it.
s/>(.+?)::(.*)/>$2_$1/': the s/old/new/ is the substitution operator and will replace old with new. Here, we look for the shortest stretch of one or more characters (.+?) from the > until the first occurrence of ::, and capture it as $1 (because the pattern is in parentheses). Everything after the first :: is then captured as $2 ((.*)). Then, we replace everything we matched with a >, then the 2nd captured pattern, an underscore and the first captured pattern.

Just for fun, here's a golfed awk approach:
$ awk -F'::' '$0=sub(/>/,"")?">"$2"_"$1:$0' file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

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


text:
PCA plot did not work in single cell RNA-seq
I run plotPCA for single cell RNA-seq data, while I get error message (I use R 4.0). I attached the code and error message here. Did anyone know the reason and how to deal with it. Thanks!
library(pcaMethods)
library(SC3)
library(scater)
library(SingleCellExperiment)
library(pheatmap)
library(mclust)
set.seed(1234567)
deng <- readRDS("deng-reads.rds")
plotPCA(deng, colour_by = "cell_type2")

Error in value[[3L]](cond) : invalid subscript 'type' in 'reducedDimNames(<SingleCellExperiment>)'
'PCA' not in 'reducedDimNames(<SingleCellExperiment>)'

I also get error message when I run the code as below:
colData(deng)$tSNE_kmeans <- as.character(kmeans(deng@reducedDims$TSNE, centers = 8)$clust)

# Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

I solve the first error for PCA, while have no idea about the second one. Does anyone have any suggestion?

The error is pretty obvious. There is no PCA data in the respective slot of the SingleCellExperiment that plotPCA could use, at least it is not called PCA which by default is expected. Use runPCA first, please read the manuals of scater and scran first, they have step-by-step code.

Thank you for your reminding. I follow a tutorial from Sanger institute and they use plotPCA(deng, colour_by = "cell_type2") without error message. While the newest scater package may have litter difference. However now I make it work and please see below: deng <- runPCA(deng)
plotReducedDim(deng, "PCA", colour_by="cell_type2")

If you follow tutorials and refer to them then please provide a link.

Here is the link: https://scrnaseq-course.cog.sanger.ac.uk/website/biological-analysis.html#clust-methods


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

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

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

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

dev.off()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Let us continue this discussion in chat.

Another 2 tree compare tools: ggtree and blatic 3

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

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

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

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

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

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


text:
Very low probability densities
I'm generating a density plot using the ggplot2 library:
ggplot(df_chr20, aes(x=start1)) + geom_density() + theme_bw()
ggplot(df_chr21, aes(x=start1)) + geom_density() + theme_bw()

The data frame for df_chr20 has 1573 rows, and I wanted to plot the probability density of genomic regions as x = start1. The same for df_chr21, with 1271 rows. What I see in the plot is what I expect, but is there a reason why the y-axis probablilty densities are so low? The area under the curve is surely not equal to 1.
I expected the high localised high density in chr 21, and both the plots have similar scales in terms of y-axis values. But since they are so low, how can I interpret this plot?
I found some information here, but my values are much lower than the examples.
Any help would be very kind!

It is not like your conventional frequency which adds up to one. The density is low because the width of your bins is huge, and the number of observations you have is low. From this by Wickham, the basic kernel is:

where K is the kernel and h is the bandwidth (Scott, 1992b)
If you don't specify it, by default h will revert to nrd0:
library(rtracklayer)
gr_narrowPeak <- import("https://www.encodeproject.org/files/ENCFF072FBH/@@download/ENCFF072FBH.bed.gz",format = "BED",extraCols = extraCols_narrowPeak)
gr_chr21 = gr_narrowPeak[seqnames(gr_narrowPeak)=="chr21"]

da = as.data.frame(gr_chr21[rep(1:length(gr_chr21),gr_chr21$signalValue)])
dens = density(da$start)
dens$bw
[1] 1133889

plot(dens)

ggplot(da,aes(x=start))+geom_density()

You can decrease h such that you get something smaller (but at the cost of looking coarse). Important thing is it doesn't sum up to one.
If you would like to explain it easier like a histogram, consider converting it into counts, like explained in this post:
binW = 50000

ggplot(da,aes(x=start))+
geom_density(aes(y = ..density..*(nrow(da)*binW)),bw=50000) +
geom_histogram(binwidth=binW,fill="steelblue",alpha=0.5) +
theme_bw()

ggplot(da,aes(x=start))+
geom_histogram(aes(y=..count../nrow(da)),binwidth=50000) +
geom_density(aes(y = ..density..*binW),bw=binW,fill="steelblue",alpha=0.5) +
theme_bw()


text:
Converting hmmer --tblout to pandas dataframe
Is there a way to convert a hmmer output to a pandas dataframe? I am also unsure how to load a hmmer tblout table into python via the Bio module. I believe you can call a hmmer format with SeqIO.parse or SeqIO.search. How can you convert the table to a pandas dataframe?

This question was answered on Stack Overflow: https://stackoverflow.com/questions/62012615/convert-a-hmmer-tblout-output-to-a-pandas-dataframe/62021471#62021471

It is rude to cross-post on multiple forums.

The manual says fields in a tabular output file are space-delimited. It should be fairly straightforward to load a space-delimited file in pandas.
From the manual:

--tblout <f> Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.

While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference.  Link-only answers can become invalid if the linked page changes. - From Review

Thank you for the pointer, @Michael. I've added an excerpt from the manual.

Thanks @RamRS ... much better :-)

This question was answered on StackOverflow:  https://stackoverflow.com/questions/62012615/convert-a-hmmer-tblout-output-to-a-pandas-dataframe/62021471#62021471


text:
Plink 1.9 equivalent of --score sum in Plink 2.0
I have previously performing Allelic scoring in Plink 1.9 using:
plink \
--score sum sumstats.txt\

However, it seems this has changed in Plink 2.0. What is the equivalent function for Allelic scoring?

plink2 --score cols=+scoresums

See https://www.cog-genomics.org/plink/2.0/general_usage#colset and https://www.cog-genomics.org/plink/2.0/score for details.

Hi a bit of explanation would be helpful

I don't have enough reputation for a comment but the place where the option to add scoresums is listed is here https://www.cog-genomics.org/plink/2.0/formats#sscore


text:
snakemake - using output file not explicitly declared as input for subsequent rule
I am writing a snakemake that, up to now should

use raw sequencing reads to produce a draft assembly using canu
use minimap2 to align the raw reads against the draft assembly, generated by canu.

The following two rules are described below:
rule canu:
    input:
        fqdir + "{sample}.fastq"
    output:
        directory("{sample}/canu")
    shell:
        "/home/bhinckel/downloads/canu-1.9/Linux-amd64/bin/canu genomeSize=4.1m -p {wildcards.sample} -d {output} -nanopore-raw {input}"

rule minimap:
    input:
        reads = fqdir + "{sample}.fastq",
        draft_assembly = "{sample}/canu/{sample}.contigs.fasta"
    output:
        "{sample}/aligned/{sample}.sam"
    conda:
        "envs/general.yaml"
    shell:
        "minimap2 -ax map-ont {input.draft_assembly} {input.reads} > {output}"

However, if I do snakemake -np I get Missing input files for rule minimap. Which does make sense, after all I am not specifying the input for rule minimap ({sample}/canu/{sample}.contigs.fasta) anywhere i.e. snakemake cannot know that the file {sample}/canu/{sample}.contigs.fasta will be generated by canu.
The thing is that canu produces several files as output (amongst which the file {sample}/canu/{sample}.contigs.fasta is one of them), that's why I specified its output as a directory.
One way of solving this I guess would be to specify all output files produced by canu as an individual rule output (e.g. {sample}/canu/{sample}.contigs.fasta, {sample}/canu/{sample}.report, {sample}/canu/{sample}.correctedReads.fasta.gz), though I am looking for a more robust way to do so (if that's possible at all).

Another way is to add one command after minimap: touch done.txt. This creates an empty file after minimap is finished. This can then be used in output rule for snakemake.
But specifying all the output files may not be a bad idea if these outputs are used in subsequent rules.


text:
How to change sequence format in my alignment file?
I have fasta file with alingned several sequences (from MUSCLE), when I open it (e.g. notepad++) they look like:

And I want dot format of identities like  and then save it in txt file.
I failed using MEGA, BioEdit, Jalview and CodonCode, how can I do this?

Hi ... this is definitely in MEGA. Could you check again? Try updating to MEGA X?

@Michael yes I tried. Export to xls and xlsx - crashed. Csv - fine, but only first 1000nt. Unfortunately I have about 40 sequences ~25kb length each. MEGA can't handle this amount of data.

okay, it must be done because it is very easy to code e.g. via AlignIO in Biopython ... you'll need to adapt it, but all you are doing is splicing through an alignment.


text:
Is there a FASTA Sequence Description Nomenclature to Indicate a Split Sequence or Subsequence?
I have a nucleotide sequence where I am interested in aligning reads distinctly to one or another portion of the sequence. However, I would like to use a visualization tool such as IGV to view these alignments in the context of the parent sequence. For example, I would like to split something along the lines of:
>seq1 parent_sequence
ACTGACTGACTGACTGACTGGTCAGTCAGTCAGTCAGTCAACACACACACACACACACAC

Into:
>seq1_1:20 parent_sequence subsequence_1_to_20
ACTGACTGACTGACTGACTG
>seq1_21:40 parent_sequence subsequence_21_to_40
GTCAGTCAGTCAGTCAGTCA
>seq1_41:60 parent_sequence subsequence_41_to_60
ACACACACACACACACACAC

Then perform alignments using the split subsequences as a reference for alignment of reads, but visualize these read alignments on the parent sequence.
While I'm sure it would be possible to write a script to manipulate all the alignments to all of the "split" subsequences back into the context of their original parent sequence, I am trying to avoid this if possible.

Is there a standard FASTA sequence description nomenclature for subsequences to accomplish what I'm describing that is recognized by programs like IGV?
If there is not, I will likely still be splitting these sequences as described. Is there a suggested standard/common method for annotating/identifying these split sequence identifiers that maximizes readability?

Edit:
The biological context for this question is analysis of RNA-Seq data collected after pulldown of a specific RNA-binding protein and digestion of overhanging read ends. The data should be enriched for a specific subsequence of interest, but will also contain the "parent" transcript. Alignments to both portions are of interest to me, but in different ways, and the alignments to the subsequence and the remainder of the parent sequence should be considered distinct. In terms of visualization it is helpful to view the alignments to all pieces in context of the "parent" transcript.
I realize that splitting the parent sequence into two subsequences for alignment will result in loss of reads straddling the two portions, but that is acceptable given the necessity of having the two sections be distinct for alignment purposes.

Dear @Dan absolutely, yes, however you would need to explain what you are working on and what your biological objectives are.

@Michael, Thanks for the comment. I am not sure why that is specifically relevant here, but I am happy to do so.

Dear @Dan, thanks biological aims and organism/field are frequently requested particularly where there is very little technical information, particularly code. It is quite common to encounter situations where the OP says, "thats not what I meant" following an answer.

@Michael That makes sense, thanks for the clarification!

IMHO there is a widely accepted standard for split reference fasta headers, there could be something IGV specific though. When I do subsequences, I do >REF_from_to.

Any further feedback @Michael?

...  in my area mixed infections are a big issue, hence my question. So if you use Illumina short-reads its complicated. If you are dealing with mRNA  .... follow @KamilSJaron advice.


text:
Adding external annotation to the phylogeny tree
I was looking to overlay annotation on my phylogeny made,which was solved here question asked where i was looking to compare two different phylogeny which i could do.
Now i want to label the tips into various groups such as if its human or canine etc.One of the possible solution i found is this one not sure how to implement it.
alignment  file
annotation file
library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

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

i want to label according to the start column which is my source in my annotation file
Any suggestion or help would be really appreciated
My working solution using this , this works as I did the alignment of sequences using full header name which kind of help to look for the pattern and grep it.
The data file
I would still like to see the ggtree way instead of putting lots of label in the plot I would only like to see if I can annotate all the tips which have one organism source as one color rather putting all the labels
pars.a
# "Phylogenetic tree with 516 tips and 514 internal nodes"

# by using the class() function
class(pars.a)
# "phylo"

# or by using the str() structure function
str(tree)
str(pars.a)
# "List of 4"
# this list includes $edge, $Nnode, $ tip.label and $edge.length
# the tree$tip.label includes family designation
#tree$tip.label  # 516 of these
pars.a$tip.label
# from the Science paper, we have seven kinase families:
# kinase categories... TK, TKL, STE, CK1, AGC, CAMK, CMGC
# with the following colours
# "red", "green", "paleblue", "orange", "yellow", "purple", "pink", "green"

# by using the grep()function on the tree$tip.label part of the object
# we can find the tip labels that include "TK/" - i.e. tyrosine kinases
#grep("TK/", tree$tip.label)  # gives a list of numbers with "TK/" in tip label
#length(grep("TK/", tree$tip.label))

grep("Human",pars.a$tip.label)
length(grep("Human",ml.a.rooted$tip.label)
)
# thus there are 94 tip labels with that are designated TK (not TKL tyrosine kinase like)

# make a vector for each tip.label called tipcol with black on all of these...
#tipcol <- rep('black', length(tree$tip.label))
tipcol <- rep('black',length(pars.a$tip.label))

# make a vector with our list of kinase categories
#kinaseCats <- c("TK/", "TKL", "STE", "CK1", "AGC", "CAMK", "CMGC", "RGC")
kinaseCats <- c("Human","Rat","Murine","Bovine","Canine","Feline","Porcine","Gallus","Bat","Infectious","Transmissible")

# make a vector of color we want:
colorsList <-c("red", "darkolivegreen3", "blue", "orange", "blueviolet", "purple", "pink", "green","brown","navyblue","purple4")

# replace colours where grep gives "TK" as red, etc in a loop
#for(i in 1:length(kinaseCats)){
 # tipcol[grep(kinaseCats[i], tree$tip.label)] <- colorsList[i]
#}

for(i in 1:length(kinaseCats)){
  tipcol[grep(kinaseCats[i], pars.a$tip.label)] <- colorsList[i]
}

# plot with edge length false to see nodes better
#plot(tree,
 #    use.edge.length = FALSE,
  #   tip.color=tipcol,
   #  cex = 0.25)

plot(pars.a,
     use.edge.length = FALSE,
     tip.color=tipcol,
     cex = 0.5)
# slow to draw due to text - a bit annoying!
nodelabels(cex=0.4)

figure

there are nice features in R ggtree package. You may use the geom_cladelabel() function for your purposes. But I haven't find yet a way to automatize clade annotation. To do it, you have to be sure of the clade homogeneity before removing tiplabs. here are some commands to do clade label annotation:
annot<-read.delim("annota.txt", header = T)

ggtr<-as.treedata(pars.a)
g<-ggtree(ggtr)
#find the node number corresponding to the common ancestor of the top clade in you tree
MRCA(ggtr, "AAL40400", "AAF25519")
#here, check that all accessions belong to the same organism/group
subset(annot, label=="AAL40400" | label=="AAF25519")
#in this case, sequences are related to bovine.
g<-g + geom_cladelabel(node=89, label="Bovine", align=F, fontsize=4)
plot(g)

#another example
MRCA(ggtr, "AAA46456", "AAF97738")
subset(annot, label=="AAA46456" | label=="AAF97738")
g<-g+ geom_cladelabel(node=73, label="Rodent", align=F, fontsize=4)
plot(g)

looks bit more complicated due to this "have to be sure of the clade homogeneity"  if there is a mix up then its difficult

you may keep tip labels where clade are heterogeneous. Or you should consider smaller clades or consider larger taxonomic groups (eg I summarized mice and rats as rodents in the second example). In any case, i think it deserve a visual/manual inspection to ensure the correctness.

okay so the better idea is to combine the larger taxa group.So first i have to find the node number then combine. Is it possible to give separate color to the clade? such as bovine one color then rodent other etc ..?okay the color can be done i saw the manual

yes it is fully customizable. Finding the appropriate number is the more fastidious step I think. There is probably a way to automatize it however.

nodelabels(cex=0.4) this is how i can label the node in base R which i was doing it is the same for the object which is made using this "g<-ggtree(ggtr)"

this a good thing to identify interesting node numbers!

Let us continue this discussion in chat.


text:
extracted counts matrix from Tximeta generalized experiment with no linked transcriptome are normalized?
I used the Tximeta vignette to construct a summarised experiment from quant files (salmon, index prepared with NCBI transcriptome).
1)Tximeta return me a non-ranged se (which I am not sure about the real implications)
because of this...
2) I used makeLinkedTxome to link the transcriptome and the gff (local files used in the index creation).
3)run tximeta on the colData to create the summarised experiment se
4) seqinfo(se) returns me a NA populated table
5) I summarise to genes (gse)
6) I extracted the gene counts counts <- assay(gse)
My questions here are;
a) Are these counts normalised?
b) If not, how can I normalise them? TMM for example? (CalcNormFactors will not work here right?)
Is there any better way to obtain a normalised count matrix from tximeta? (maybe I don't need to link the transcriptome...)
Any help would be appreciated!
Thank you! and regards,

tximeta will summarize transcript abundance estimates (or counts) to the gene level and produce raw counts. All these metadata can be used and might be useful in certain situations, but none of these are essential for the analysis. You can equivalently just use tximport to produce the same raw counts. No, counts are not normalized. TMM should work just fine.

Thank you! so I can use these `counts'  with the pipeline. Perfect! thanks!!!


text:
Error: Insufficient values in manual scale. 24 needed but only 1 provided?
folks.
Can you give me idea about an issue that I encountered?
I ran single-cell RNA sequencing data using Suerat in R, and when I tried to draw violin plot, there is an error as below.

VlnPlot(tumors, "Cd79a", "Foxp3", "CD14")

Error: Insufficient values in manual scale. 24 needed but only 1 provided.
Run rlang::last_error() to see where the error occurred.
However, I can draw the violin plot for Cd79a. The problem happened to only several genes. I haven't solved this problem yet. Do you know what should we do?
Thank you.

You must supply genes (or features) as a vector, not individually:
VlnPlot(tumors, c("Cd79a", "Foxp3", "CD14"))

By the way it is probably Cd14, not CD14.

Thank you. I learned a lot from you. Such a stupid mistake:)


text:
GSEA preranked different results opposing phenotypes
I am running an analysis of two phenotypes: mutant and control.
I have a statistic used to rank the different genes, and the ranking is exactly the opposite between mutant and control. Moreover, there are no duplicates in the ranking (which excludes that genes are ranked alphabetically). I uploaded the test file here.
My expectation is that running GSEA using mutant as positive phenotype, should give me the same results (but in opposite direction) than running GSEA using control as a positive phenotype - since there are only two phenotypes, and the ranking statistic is computed mutant vs control.
However, when I run Gene Set Enrichment Analysis using mutant as positive phenotype, the results I get are not the same with respect to running GSEA using mutant as negative phenotype.
I initially thought that it might be related to the choice of the random seed, however, I get the same results even setting the same random seed in all runs.
The following is an example of an enriched result of the positive phenotype:

The same results with the negative phenotype:

The ranking is exactly the same, however, the leading edge is different in the two cases. It's probably due to the fact that the ranking sum goes in one direction (increasing the enrichment score) in both phenotypes. This is probably the cause of different scores in the opposite phenotypes, however, which is the result I should use in this case?

@ATpoint I have added an example of the results because I think it is useful to identify the issue here. I also added the ranking test file. My code is just the execution of the Broad version of GSEA, I can provide the data if necessary, but I think that the results example I added should be sufficient. Does it work?

There are two things that I notice from these plots that are unusual for GSEA analyses that people do:

The zero cross is at the extreme end of the plot, meaning that all selected genes have a positive test statistic. I think it's a good idea to separate out the positive and negative test statistics (because then the zero point cross is at zero on the Y axis), but other people disagree with me, preferring that the plots show all genes, rather than positive or negative genes. In any case, GSEA seems to work best when the ranked gene list is unfiltered.

The "negative phenotype" plot has an unusual shape, with high enrichment score boosts for what should be low-ranked genes (and, now that I look at it, the first "positive phenotype" plot has similar issues). Most often, I see GSEA plots that have gradually decreasing enrichment score boosts as the graph moves further down the list of genes. The "negative phenotype" plot also doesn't have any sign changes; the bottom bits should be blue, not red, if the flip has been done correctly.

I suspect what's happened here is something like that there is no test statistic in the list of ranked genes, and the gene rank is being used as a distinguishing statistic. For the GSEA to look "normal", the gene list should have a field that indicates signed difference (e.g. mean expression of mutant minus mean expression of control). With a list presented in this fashion with signed difference, it doesn't really make sense to compare mutant and control phenotypes separately, as the test statistic incorporates both at once.

What do you mean by having the same results but opposite direction?
What GSEA does is that it goes through the gene list from the top to bottom, whenever it encounters a gene that belongs to a gene set, it adds a positive score to that gene set, and if not, it penalizes the score. That is what you see in the green line in the results.
Now imagine that you reverse the order of your input gene list, it is obvious the results are going to be different if the program still goes through the list from top to bottom.
I think as a general rule, whenever you change the order of input gene list, results are going to be different.

There's a misunderstanding. With 2 phenotypes, e.g. Case and Control, the ranking is the opposite, correct. However, I can test Case by looking at genes upregulated using the "Case" ranking (top genes upregulated in Case), or by looking at genes downregulated in Control. GSEA provides the two results for each phenotype (e.g. genes upregulated and downregulated in Case). I understand what the running sum does, and I think the difference in results (for each phenotype, e.g. Case) is due to the nature of the running sum and the leading edge. What I don't understand is which result to consider.


text:
How can I plot multiple sequence alignment in R as a heatmap
I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.
I came across to msaplot from ggtree package. I also checked ggmsa package for that. But so far, I did not get what I wanted. So basically I wanted to:

change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)
plot as a heatmap with two colors based if the amino acid changed or not

here's some data example of what I doing, just instal ggmsa and you can try this:
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)

Could you post your code?

I added a very simple example code present on ggmsa package

@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.
@StupidWolf:
We read in the alignment:
library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300)

Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :
aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:
bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:
library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.


text:
Converting aligned fasta to plink ped/bed
I have an alignment of multiple sequences in a FASTA file (output from MAFFT), for which I would like to simulate a phenotype using plink, but for that I need to have my alignment in a PED file or plink BED. There's plenty online about converting in the other direction, but the only mention I found of converting FASTA to BED is about aligning a raw FASTQ file for use in plink - my data is already aligned! How can I do this conversion?
Note that How to convert FASTA to BED is about the UCSC BED format, not the plink binary PED.

FASTA file contains sequences, whereas PED file contains variants. Usually, you need to do variant-calling on sequencing file to get variant file. I do not know that is contained in your FASTA file.

I have the same challenge, I have aligned sequences but not vcfs, bams or fastqs. Any ideas?

An alignment can be the result of two slightly different analyses:
There is multiple sequence alignment (which is what you get from MAFFT) where sequences are aligned so that similar regions are on top of each other. This may require introducing indels (insertions / deletions) if a particular region is absent in some of the sequences.
Then there is alignment to a reference sequence. This is usually how you get variant data and most of the tools deal with this use case. Therefore, in the link you provided the solution includes aligning the data to a reference sequence. Usually the refernce sequence is much larger (e.g. the human genome) and after alignment the result is a file (.bam) that tells you where the query sequences match in the reference.
I did a quick search for ways of converting a multiple sequence alignment (MSA) to VCF. There is a tool called msa2vcf in the Jvarkit collection of utilities that can do this. The example is for the CLUSTAW format but FASTA is accepted as well.
You shoud check that indels were treated correctly because these are most likely to cause trouble.
Then you can simply convert the VCF to PLINK format  using PLINK.


text:
Clustering issue about overlapping CD4, CD8, and Foxp3 genes in mouse tumor models
bioinformatics scientists.
I hope your work goes smoothly well and be safe.
I'm stuck in trouble at the step of clustering.
Briefly, I downloaded fastq files from a published study (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861). There are 24 samples that come from 6 different mouse tumors.
Each tumor consists of two biological replicates where one biological replicate includes two technical replicates. So, I ran the command "cellranger count", and generated 24 "filtered_feature_bc_matrix" files. For example, for the MC38, bio1-tech1, bio1-tech2, bio2-tech1, bio2-tech2 samples exist. FYI, These samples contain tumor cells as well as tumor-infiltrating immune cells.
And then, I carried out the command "cellranger aggr" to combine all matrix files into one matrix file.
As suggested Seurat tutorial (https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html), Standard pre-processing workflow, linear dimensional reduction, determining dimensionality, clustering, and differential expression were done. Below are the commands I put.

tumors <- NormalizeData(tumors)
tumors <- FindVariableFeatures(tumors, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(tumors)
tumors <- ScaleData(tumors, features = all.genes)
tumors <- RunPCA(tumors, features = VariableFeatures(object = tumors))
tumors <- JackStraw(tumors, num.replicate = 100)
tumors <- ScoreJackStraw(tumors, dims = 1:20)
ElbowPlot(tumors)
tumors <- FindNeighbors(tumors, dims = 1:10) # I chose dims from 5, 6, 7, 8, 9, 10 or 20 as a cutoff.
tumors <- FindClusters(tumors, resolution = 1.2)
tumors <- RunUMAP(tumors, dims = 1:10)
VlnPlot(tumors, features = c("Cd79a", "Foxp3", "Cd3e", "Cd4", "Zap70", "Cd8a", "Cd14", "Il2ra"))

Interestingly, as you can see here, Foxp3, Cd4, Cd8a are found to be overlapped. This phenomenon was the same when I drew vlnplot using different dimensionality from dims=1:5 to 1:10 or 1:20.
Next, I performed this single-cell RNA sequencing analysis using each tumor; mc38, b16, sa1n, ll2, emt6, and ct26. Representatively, MC38, SA1N, EMT6, and CT26 are presented as below.
For mc38, CD4, CD8, and Foxp3 are overlapped.

For sa1n, CD4, CD8, and Foxp3 are separated.

For emt6, CD4, CD8, and Foxp3 are overlapped.

For ct26, CD4, CD8, and Foxp3 are overlapped.

So, except for SA1N that only separates Foxp3 from CD8 T cells, I found that other tumor types can't separate Fopx3 from CD8 T cells. Do you think it is related to technical issues, such as values of dimensionality or PCA I put? Since I am not an expert about scRNA, I don't know the exact meaning of dims or PCA, though. Is there anyone who can tell me what's going on? And, do you know how to make the right clustering to get exact T cell subpopulations?
I do appreciate your valuable time and thank you for reading this long story.

It could have something to do with the parameter you chose for clustering. Usually, people try many sets of parameters. Also, it is somewhat known that CD4,CD8, and Treg can be really hard to be separated.

You can check where the CD4, CD8, and Foxp3 clusters are. The short explanation is that the genes that you are interested in, are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs

Then it's a matter of deciding the resolution to find the clusters.. You can try and play around with this.. I would only have time tomorrow to look at this data and write a solution

@StupidWolf, Hello again. I solved this issue by increasing the resolution up to 1.8 or even higher. I also happened to change the value of min.cell or min.features or others, such as nFeature_RNA, nFeature_RNA, but it appeared that these are "not that critical" enough to change the result. Your advice sounds really important, but I don't understand the part "are weighed differently in the PCs, so you need to first checked that out, and be sure to include those PCs". Are PCs the PCA? Do you mind If I have a basic question about the terms?  What is the difference between PCA and dimensionality?

You can check this for a visual guide to PCs, https://setosa.io/ev/principal-component-analysis/

hi @raiora, glad you manage to find the cluster ! sorry i haven't found time yet to look at it.. Yes the PCs are the results from PCA. You have a lot of genes, and you can do what is called dimension reduction as you have pointed to. PCA is a method of dimension reduction which finds the linear combination of genes that explains the most variance in your data. The data is projected onto so called principal components which is what I referred to as PCs. This is what you are looking at with the elbow plot, how much variance each PCs explain

So from the first step, RunPCA(tumors, features = VariableFeatures(object = tumors)), you can check whether foxp3 etc is in VariableFeatures(object = tumors)), because this will decide whether the PCs capture the variation in t-reg cells. By default i think 2000 genes are used? Then you can look under your seurat object, tumours@reductions$pca@feature.loadings to see what is the loading for foxp3 etc


text:
MEME's motifs have an E-vaue too high
I'm currently learning how to use MEME-Suite (web tool) to find a specific known transcriptor binding site in the promoter region of 5 putative bacterial genes of interest, but the results obtained with just the five promtoer sequences have an E-value too high. What can I do to improve these results? Can I use as input some other reported sequences that have the motif I'm looking for without introducing any bias?
I used the option ZOOPS, with a background model of first order, and I limited the search to palindromes only, since I'm looking specifically for a 20 bp sequence with dyad symmetry. The input sequences are the promoter regions of each gene (200 pb upstream of the start codon).


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:
GWAS Rooted PCA analysis problem
I'm fairly new to plink software and wanted to get some additional practice after doing several tutorials. I obtained the data from this paper (I'm not using this paper's methods) to do some QC with subsequent analysis because they mostly use SNPstats.
Some background on what I did leading up to PCA:
My QC steps involved:
Filtering for low MAF
Filtering for high missing values on the individual and SNP level
Excluding outliers in heterozygosity (based on F statistic)
Excluding controls that deviated from Hardy-Weinberg
Finding related individuals and excluding members from the pairs that had the lower call rate
After these steps, I had 657649 variants remaining, and no individuals had been filtered out from the original total (1401)
I then used a 1000 genomes dataset that I had already downloaded (unsure of which phase, but it contained 690 individuals) and merged it with my data (did similar QC steps on the 1000 genomes data, and also used the --flip flag to fix strand issues).
I did a rooted PCA analysis (with a LD threshold of 0.2, which included 112795 variants post-filter, but I also tried PCA without the filter with similar results) to look at population substructure and got this plot (in R):
Where the pink "OWN' labels in the plot are my data points. Does anyone have any idea why my data wouldn't be lining up with any of the ethnic groups from 1000 genomes? The data is said by the authors to be "pre-filtered," which I'm assuming means fairly homogenous, so I wasn't expecting them to be in different ethnic groups, but they should be lining up with one of them, right?

It would be helpful if you attach your resulting plot here.

Agreed 100% with @PhoemixMu


text:
about the scaling after normalization and transformation in RNA-seq data
When we use count data in RNA-seq analysis, we usually use normalization and sometimes vst, rlog transformation (DESeq2)or log2(CMP+4) transformation (edgeR) to perform K-means clustering. Can we use scaling after normalization and transformation?
Or Should we use scaling only when we just use normalization without transformation?


text:
Motivation behind the neighbor-joining distance matrix recomputation
In each iteration of the neighbour-joining method of phylogenetic trees construction, after joining the nearest neighbours the (additive) distance matrix $D_{m \times m}$ is recomputed as
$$Q_{ij} = (m-2)d_{ij} - \Bigg( \sum_{k \neq i} d_{ik} \Bigg) - \Bigg( \sum_{k \neq j} d_{jk} \Bigg)$$
As is stated in numerous sourced that I found, it prevents the joining of non-direct neighbours, but am looking for a more detailed explanation and motivation.
Sources: Wikipedia

Let's look at the components of the equation one by one first since it is easier to describe each part on its own. We can do this because they are added together in computing $Q_{ij}$. Component1, $(m-2)d_{ij}$ is saying that we want a quantity which would be the total penalization, if every distance for $m$ sequences was $d_{ij}$ excluding the 2 distances between themselves (it's a direct proportionality so that the $Q_{ij} \propto (m-2)d_{ij}$). Component2,
$\left(\sum_{k\neq i}d_{ik}\right)$, says that we look at the sum of the distances sequence $i$ has against all other sequences excluding itself
(it's an inverse proportionality so that the $Q_{ij} \propto \left(\sum_{k\neq i}d_{ik}\right)^{-1}$), and $\left(\sum_{k\neq j}d_{jk}\right)$ as component3 is analogous.
non-direct neighbours means that on the phylogenetic tree, 2 extant sequences (sequences in the alignment) have another sequence set that can link them together. If two sequences are definitely related to each other then on the phylogenetic tree they will produce a bifurcation and the smallest distance between themselves as being directly related, $(m-2)d_{ij} < \left(\sum_{k\neq i}d_{ik} + \sum_{k\neq j}d_{jk}\right)$, eg. the $d_{ij}=1$ and the rest are all greater than 1 gives us a low number. If $f_{max_{i,j}}(d_{ij})$ is used we'd be using the maximum distance for component1 so the overall measure would be positive and as large as possible.
In summary it provides a measure for placing the sequences together with the closest distance. You may say that this is 'overkill' as we can just search the $m^2/2-m$ set of distances and choose the minimum. In the goal for producing the tree, things get more interesting. Let's say look at 3 sequences $i,j,k$ of a larger set $m$, and see that sequences together have the same distance between them $d_{ij}=d_{ik}$ in respect to $i$. $d_{ij}=d_{j*}$ would mean that $j$ is equally distance to all other sequences as it is to $i$ and let's say $d_{kz}=max \land d_{ky,y\neq z}=d_{ij}$ (the sequence $k$ is equidistance to all sequences except to a particular sequence $z$ and can happen in some region in the alignment). Which do we join with $i$? According to the formula we choose $k$ over $j$ because $k$ is more distant from other sequence(s) as it would provide a more negative number (lower values cause of the 2nd and 3rd components although component1 is the same in both cases). Effectively we want the minimum distance to $i$ and simultaneously the maximal distance from other sequences (penalize sequences which could be closely related to the other sequences helps produce an overall better tree).
Anders Gorm Pedersen (of DTU) in his coursera videos (also on YT) provides some numerical examples to see how the numbers are produced on an example tree.

NJ as with all clustering methods there isn't an explanation, It is a hierarchical clustering method but there are many algorithms to cluster a distance matrix. NJ is closely associated with the concept behind UPGMA (unweighted pair group mean algorithm)
The key difference with UPGMA is that after joining two taxa by NJ the reconstructed genetic distance between the node and one taxa versus the other taxa is not assumed to be equal, i.e. the rate of mutation is heterogeneous. All other clustering algorithms I know assume a "global molecular clock" if that makes any sense, the mutation is constant against time. NJ removes that assumption and thats why we like it.

So there is no motivation behind the formula other than 'it just works'?

Well, clustering is not statistically cool no how you look at it. HOWEVER, NJ allows the calculation of a additive tree  (phylogram) rather than an ultrametric tree (cladogram) and that is very useful.


text:
separate genes with expression values
I have a file like this with 2 columns. I want to separate the columns based on plus and minus signs into separate files. How do I do this with awk or sed. Kindly help.

What do you want to do with NA values?
The below scripts ignore them (h3ab74's answer will print them to positiveExpression)
awk '$2 ~ /^[0-9]+$/  {print $0}' expressionFile.txt > positiveExpression.txt
awk '$2 ~ /^-[0-9]+$/ {print $0}' expressionFile.txt > negativeExpression.txt

These do a regex match to respectively, positive and negative numbers. If there's anything else in field 2 the line won't be printed.

awk '{if($2>0) print $0}' expressionFile.txt >> positiveExpression.txt
awk '{if($2<0) print $0}' expressionFile.txt >> negativeExpression.txt

Assuming that the file with gene expressions has the title: expressionFile.txt, this should do the trick.
Explanation of first line of code: if values in the 2nd column - $2 - are greater than 0, print all columns ($0) and write to a file (>>) titled, positiveExpression.txt. The use of awk '{}' is just syntactical.


text:
apply TMM on counts imported from salmon using Tximeta
I used Tximeta to import a summarisedExperiment from the salmon output (used with genocide transcriptome v34). I need to produce 4 matrix of counts:
- tx counts in TPM
- gene counts in TPM
- tx counts normalised
- gene counts normalised
The first two I used assay(se, "counts") and assay(gse, "counts") respectively.
Mu problem is to apply the normalisation to them. I used deseq2 to estimate the size factors with a fake design because there are no groups yet in the dataset so I used ~gender to been able to generate a adds object to apply the EstimateSizeFactor, but when I try to see them returns me 'NULL'
> dds_genes <-DESeqDataSet(gse, design = ~ gender)
> dss_genes <- DESeq2::estimateSizeFactors(dss_genes)

> DESeq2::sizeFactors(dss_genes)
[1] NULL

I get count values like this after using these commands (which I am not really sure they are really normalised;
normalized_counts <- DESeq2::counts(dss_genes, normalized=T)
normalized_counts[1:10, 1:10]

ENSG00000000003.15  654.92903  771.34051  500.1863  560.30530 1144.8295  938.40446 1173.208233  896.28964  388.1482   501.20545
ENSG00000000005.6     0.00000    0.00000    0.0000    0.00000    0.0000    0.00000    1.344747    0.00000    0.0000     0.00000
ENSG00000000419.12  416.59305  399.05434  675.0911  285.64627  551.8253  424.59423  461.448691  679.08099  365.8343   461.66276
ENSG00000000457.14  242.52226  244.91792  279.9109  414.24494  262.8790  354.50834  451.671501  427.08499  407.7070   329.50667

I don't know how to proceed to just get a TMM normalised values from tximeta object. That would be my preferred option. I guess deseq2 can also provide me with normalised data.
Maybe there are other ways?
Thank you for the help!!
Regards


text:
about rescaled data (median)
When we have a rescaled data with median (divided by median value)(e.g.,metabolome) and perform natural log transformation, do we need to use autoscaling again to perform PCA or to compute euclidian distance?

I understand that PCA and computing Euclidian distance need scaled data. Do I need to autoscale again to these median-rescaled data (simply divided by median), Some companies give us rescaled proteome data.

Voting to close since you already asked this. Please add details to the previous question instead of spawning similar questions. https://bioinformatics.stackexchange.com/questions/13424/about-the-scaling-after-normalization-and-transformation-in-rna-seq-data

I would like to know when we have median scaled data (divided by median) in metabolomics , some use natural log transformation and to compute euclidian distance and others use natural log transformation with subsequent autoscaling to compute euclidian distance. So I would like to know which methods we should use.

Whether to scale your data again in the PCA function even if you have already had some reasonable ways to scale it beforehand has no clear answer. Theoretically, there is no clear reason to do it or not. I guess you could try both and look at the resulting figure and decide.


text:
Compare nucleotide mutation with amino acid
To begin with I have two nucleotide sequence to compare I want to see how they might have evolved in terms of sequence similarity and differences. One of the way i can think of is align the sequences. But i m not getting how do i associate it with mutation.
Now question is how do i know about mutation. I have the respective amino acid data from which i know there been mutation which is making changes how the binding happens between two sequences.
So now i have to associate nucleotide mutation with amino acid sequence data.
Any suggestion or help would be really appreciated. As always if there is a R way to performing the analysis then it would be really good.

Are you familiar with non-synonymous vs synonymous analysis and is this part of the phangorn package?

i did come to read that let me explore it..will get back

Thats probably what you want ... depends on your investigation

so my simple idea is compare the respective nucleotide sequence and see if there exist a mutation if yes then do the same in the protein sequence. So to see mutation I'm not sure how to start ..alignment perhaps not sure though..how to approach it..

You must have done an alignment to generate the phangorn tree and if you have branches and stuff you must have mutations. An alignment viewer that converts in and out of amino acid translation is a place to start

yes i have done mutation and i have the tree as well ."An alignment viewer that converts in and out of amino acid translation" jalview would be a good start i guess..

Cool, I will sketch a few ideas below when I have a minute

that would be really kind of you

I  would output the aligned nucleotide sequences and input them into an alignment editor such as Sequtron (Mac), MEGAX(particularly PC) or Jalview (Linux).
The alignment editors will allow you to flick between the triplet codon and the amino acid sequence. Therefore spotting which triplet codon mutations associate with each given amino acid mutation is easy.There is no R equivalent that I am aware of.
More generally I would do a non-synonymous / synonymous analysis which fangorn will perform. This assesses the number of mutations in which cause an amino acid change against the number of mutations that cause a 'silent' mutation (no amino acid change). A synonymous phylogeny will represent all mutations which don't change the tree. A non-synonymous phylogeny is the same as an amino tree.
There are character mapping approaches which will map an amino acid mutation onto a phylogeny as an alternative.

I did run the kaks but ran into this https://bioinformatics.stackexchange.com/questions/13499/calculating-synonymous-and-non-synonymous

Hi @krushnachChandra I've answered your question. Why it fails is easy, but correcting it is more complication, however this is needed before any form of amino acid to nucleotide mutation analysis can take place.


text:
List of global ethnicities
Looking at the 1000 genomes project I see the following list of population groups and super-groups:
https://www.internationalgenome.org/category/population/
However, is there something like a standard list of global ethnicity? Like an ontology of population diversity or similar?
I'm trying to tie genome data to 'ancestry', and to do that I'd like to survey a set of people on their ancestral background (as far back as grandparents ideally).
How to capture diversity in terms of what people 'self report' in their backgrounds? e.g. middle / far east, north / south India, etc.

It's not clear exactly what your goal is. It might help to have an idea of how you are trying to use this data; and it's probably not a bioinformatics question, as ethnicities are social/political rather than biological constructs.
However, there exist some controlled vocabularies used by LoC for ethnic groups; that looks a bit old so it might not be the best source.
I am by no means an anthropologist but it does look like there are anthropological databases devoted to e.g. ethnography that may have complementary information.
Wikipedia has a simple list of ethnicities that notes, instructively, "There has been constant debate over the classification of ethnic groups."
I think that you are unlikely to find a single authoritative answer that is safe to use as a source of truth. I would normally suggest to try another stackexchange site, though it looks like there is none for sociology or anthropology. There is no anthropology tag available on this site, but you might try posting in e.g. history.stackexchange.com and adding an anthro tag.

Thanks, tried to improve the question a bit. Some good links :-)


text:
How to calculate mutation rate and mutation sites in a genome using FASTA file?
I have 6 viral genome sequences of the same virus and 1 reference sequence in FASTA format.

How I can identify mutations and mutation sites in those genomes using FASTA sequences  but how I can do this for FASTA file?
And how I can identify the mutation rate for one genome?

Background If I've FASTQ file then I'll simply align the reads to the reference and by using variant calling tool I will get the mutate sites

The mutation rate is impossible to calculate without knowing when the samples diverged. If you have fasta sequences, you can align them and build a tree (I guess beast is the most common tool for viral genomes including time calibration).

Typically you can align FASTA or FASTQ sequences.  If your variant caller of choice doesn't accept a FASTA and you trust the bases in your FASTA you could mock up the quality scores and make a FASTQ

You can still align a FASTA file with a tool like bwa mem ( if they are short reads ) or minimap2 for long reads and run it through a variant caller like freebayes. Alternatively, if you have sufficient reads to create an assembly of all of the genomes, or if they already are you can create a MSA of them using a tool like Cactus (as it is built for genomes) or MAFFT if the sequences are short enough and then use your favorite viewer to find variable sites

I've assembled genome and the genome size is about 30-32k, so I align them and after that, I try to check variable region, It will take a lot of time.

And how I can calculate mutation rates?

I don't expect it to take too long to run as MAFFT i-INS-n ( mush slower then Cactus but probably more accurate)  can align my dateset of 208 seq at 200kb  in around 20000 clock hours  (20 hrs with 10 cores ) so 6 should be fine unless they are really gappy.  I am unclear what you mean about mutation rates in this context unless these genomes have established time points


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.

I personally rank by signed nominal p-value. The fgsea author recommended on biostars.org to use either -log10(nominal p-value) or the F-statistics column (or whatever statistic the tool you use outputs) but not FDR/adjusted p-values as the latter produces a lot of ties for genes with low significances.

Please note that I edited the post. It must of course be -log10(nominal P), no nominal P itself, sorry.

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


text:
mining DNA barcodes from Genebank/BOLD per location
Is there a method to download all the sequences for a particular geographic region (Panama) for 'DNA barcode' sequences (ITS2, rbcL or trnL). Hopefully, the specimen collection location will include this information.
I want to use the Genbank sequences to create a DNA reference database of local plants for eDNA studies.

Hi @uvnomad yes its easy, but what does your reference database do? Is it for taxonomy of all plants?

Great, how is it done? The reference is used to compare sequences of unknown plants (eDNA) to the known sequences of locally occuring plants.


text:
What are the p01sum and p0prop parameters in godon?
When using branch-site model in Godon, model optimization results include the following parameters: p01sum and p0prop. What are those and how do they relate to p0 and p1 in the branch-site model?

Disclaimer: I'm the author of godon. The question above comes from the real communication with a software user.
In the godon implementation of the branch-site model I performed the following reparametrization:
$p_{01sum} = p_0 + p_1$
$p_{0prop} = \frac{p_0}{p_0 + p_1}$.
This has the following advantages:

Both parameters have a well defined range $(0, 1)$.
There is less dependency between the two. E.g., $p_2=1-p_0-p_1$ (proportion of sites under positivie selection) does not depend on $p_{0prop}$.
Both properties are very helpful for both likelihood maximization and MCMC.

In case you are interested in $p_0$ and $p_1$ there is a straightforward way to go back to the original parameters:
$p_0=p_{0prop}*p_{01sum}$
$p_1=p_{01sum}-p_0$.
P.S. I also added this information to the tutorial.

Damn, I totally missed this package. Godon is seriously nice, great work!


text:
Determine if an organism is microbic or unicellular
I am doing some metagenomic analysis of samples which can only contain microbes due to the experimental setup.
For this, it would be useful to be able to distinguish whether a specific organism is a microbe or not, as it allows:

to detect potentially false findings,
to avoid an unnecessary inflation of the metagenomic database.

For animals and plants, this is straightforward, as they are never microbes.
However, for fungi (and some other clades) it gets complicated.
At the very least, I failed to find a simple list of microbial or non-microbial clades.
Ideally, I would like to have a function that maps a taxonomic ID to a boolean stating whether this organism is microbial. I acknowledge that this is utopic since microbe is not that well defined, there are facultatively multicellular organisms, and for some cases it may simply not be known. I am therefore happy with any reasonable approximation such as unicellularity. In particular, given my application, a small chance of falsely identifying an organism as microbial is not an issue.
I am currently working with the NCBI database accessed via the ETE 3 module for Python, but I am open to use other databases and tools.

Hi, you may eventually try to find marker genes related to intercellular communication? I do not know whether there are proteins specific to algae synapsis or cell junctions in fungi for example. Howver, it may be quite challenging to identify such markers!

@thomasdugedebernonville: That multicellularity has evolved roughly thirty times and lost a few times, I doubt that anything universal exists. Cataloguing the tree of life by marking the nodes where multicellularity emerged and ceased certainly beauts that.

Yes but convergent evolution probably drives the appearance of conserved functions at cell boundaries. However, proteins with similar roles will probably too divergent for your purposes. By extending your statement " anything universal exists", it is gonna be pretty tricky to get  a reasonable approximation for your classification. Hope you will be able too, however!

@thomasdugedebernonville: Mind that my point is not that I want something universal, but since there is no reason to expect something universal, I would have to check all cases of multicellularity anyway, in which case I might as well build a database of those myself.

makes sense. I can be a very nice input to pre-existing databases!

you may have a look at this publication maybe? https://academic.oup.com/gbe/article/8/4/1279/2574110

I would start by saying that your statement that no animals/plants are microbes is false. There are indeed microbial algae (plants) and microbial animals.
So it is indeed difficult. (I'll briefly note that unless your protocol specifically gets rid of environmental DNA, you will probably have some contaminating non-microbial DNA no matter what.)
However, we can use a different version of your proposition to make a simple (over-simple, most would argue) approximation: all prokaryotic organisms, e.g. bacteria and archaea, are microbial, and all eukaryotic organisms are not. This is not hard to do if you are already using NCBI taxonomy, as you can use e.g. this mapping. It has the advantage that all bacteria and archaea are definitely microbial.
Of course it's quite wrong, as you already know. In principle you could set up a complex decision tree sorting this out, though I imagine that's the thing you wish to avoid.
In terms of other genome databases, you might consider looking into GOLD, which has extensive metadata in a slightly less arcane setup than NCBI, which might be helpful. But there is a high rate of missingness there, and also there is still not necessarily a great "microbe"/"not microbe" mapping.
There are also existing databases of microorganisms, which may or may not satisfy your requirements; I have little familiarity with most of those so caveat emptor.
I believe (hope) that NCBI taxid is a good way to do lookups in all such databases, but it may not be true.

Thank you for your answer, but your simple approximation is indeed too simple as the vast majority of eukaryotes in my samples (and probably others with similar issue) are unicellular microbes, mostly fungi. — In principle you could set up a complex decision tree sorting this out, though I imagine that's the thing you wish to avoid. – Well, if I go by unicellularity, I would “only” need a list of the roughly thirty clades that developed multicellularity and the few cases which lost it. So far, this seems like the most simple solution to me.

@Wrzlprmft if you feel confident to implement this decision tree, then I think that might be the best approach. My concern in that case is that basal eukaryotes are relatively quite poorly described and that might affect your results, but I don't have a better alternative. If someone hasn't already implemented a procedure like this, it would honestly be a bit of a public service to the metagenomics community for such an approach to be publicly available.


text:
What is the full script to extract sequences from fasta file by using Ids in text file in python 3 and pycharm?
I have a fasta file (gene.fasta)with sequences with names patterns:
>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA
TTGCTCATTTTGAATGCCGTTACATTGGAGAAGGATGCTAATGATTATGCCGTTGGCGAT

and I have the Ids in text file (seq.txt) which are not the exact as in the fasta file:
HSC_gene_996
HSC_gene_9734
and some of the names came as
HSC_gene_996|HSC_gene_9734

How can extract the sequences?
I am not experienced in python so please use python for dummies language :)
Thank you in advance

Does it have to be python? There are far easier ways using only bash, etc.

I just have python in my computer, if there is please help me, I tried to find some online tools but it didn't work

You may want to try Joe Healey's bioinfo-tools repository. This script looks like the one to use.

To second @RamRS, that script should work indeed - provided your IDs are exact full or substring matches (excluding the >). It's been a while since I used it, so if it could be tweaked to be more useful, I'm happy to take suggestions. Whilst I love bash etc., I'd probably disagree with NatWH that it is easier than using something like biopython which is purpose built for this.

OP's problem seems to be their platform. Python excels in being a cross-platform solution, hence OP looking for Python based solutions makes sense.

Hi @SamerHabash, i can write it in python. Can you give an example where it is supposed to hit (i suppose Hsc_gene_9627.t1 and Hsc_gene_9627) and instances where it doesn't?

Hi @StupidWolf, that would be great thank you. All my target sequences IDs are like (HSC_gene_996) in a text file which are missing the (.t1 and .t2 in some cases). in the fasta file they are all ending with ( .t1 or some cases .t2- like this one :  >Hsc_gene_9627.t1 )

For example, the fasta are like these:
gene.fasta:

>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACA
ACATCTGCTGTGCAGCCGCTTCA
>Hsc_gene_9627.t2
GGGGTTTTCCCC
>Hsc_gene_962.t1
AGTCGTCAGTCAGTAGTCGC

seq.txt:

HSC_gene_9627
HSC_gene_9734
HSC_gene_996|HSC_gene_9734

We use a few modules, read in the fasta:
import re
from Bio import SeqIO
records = list(SeqIO.parse("gene.fasta","fasta"))

Read in the list:
seq = [i.rstrip().lower() for i in open('seq.txt').readlines()]

Here I use two functions, one to write fasta files, and the other to format the header for comparison:
def fasta_out(rec):
    return(">"+str(rec.description)+"\n"+str(rec.seq))

def format(rec):
    return re.sub("[.][^ ]*$","",rec.description).lower()

Define the records to keep:
keep = [fasta_out(i) for i in records if format(i) in seq]

Write them out:
f=open('matched.fasta','w')
f.writelines("\n".join(keep))
f.close()

We can check:
open('matched.fasta').readlines()
Out[50]:
['>Hsc_gene_9627.t1\n',
 'ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA\n',
 '>Hsc_gene_9627.t2\n',
 'GGGGTTTTCCCC']

Most likely there are some simple while, for loops that does not require modules..

Thank you very much @StupidWolf


text:
What steps should I follow for DNA analysis, for a classification problem?
I have a bed file which contains DNA sequences information as follow:
**
track name="194" description="194 methylation (sites)" color=0,60,120 useScore=1
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +
chr1    714176  714178  FALSE   12  -
chr1    720864  720866  FALSE   988 -

**
I loaded the bed file in R and named the matrix DataSet.
I used the follow code to get the sequences:
mydataSet_Test1<-dataSet[,1:3]

library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet_Test1[,1]))

#50
start<-as.matrix(as.integer(as.character(mydataSet_Test1[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet_Test1[,3]))+50)
Seqs50_Test1<-getSeq(genome,chr,start=start,end=end)

I'm having DNA sequences and for each sequence, I have its Methylation Number (i.e. 800). I have split the data to hypomethylated (methylation < 200, class 0) and hypermethylated (methylation > 800, class 1) and I want to build a model that can predict if a sequence (or a group of sequences) is hypo or hyper - methylated.
I now have a matrix of 358.367 data and 2 columns. The first column is the
DNA sequences and the second column is the class for each sequence (0 or 1).
I want to build a classification model in R, using SVM algorithm and 83 features (dinucleotides, trinucleotides, etc.).
Before calculating the features and building the model,
are there any steps I should follow for my DNA sequence data (each row of the matrix is a sequence from the human genome),
like filtering CpG islands or sequences, homolog reduction or overlapping sequences etc?
As you understand I have a lack in DNA analysis.

Hi Giannis, could you please describe the classification problem you are trying to solve?

I'm having DNA sequences and for each sequence, I have its Methylation Number (i.e. 800). I have split the data to hypomethylated (methylation < 200, class 0) and hypermethylated (methylation > 800, class 1) and I want to build a model that can predict if a sequence (or a group of sequences) is hypo or hyper - methylated

@GiannisLazaridis Please [edit] the question to include the last comment, and include the code or guides you have tried

I edited the question


text:
about the cpm normalization after using normalization factor
Is it okay to use CPM normalization (with/without log transform) after using TMM normalization? Why do we need both?
 library(edgeR)
 library(SummarizedExperiment)
 load(url("http://duffel.rail.bio/recount/SRP049355/rse_gene.Rdata"))
  counts <- assays(rse_gene)$counts
  y <- as.matrix((counts))
 y <- DGEList(counts = y, group=c(1,2,3,4,5,6,7,8,9,10))
  y <- calcNormFactors(y)
 z <- cpm(y, normalized.lib.size=TRUE)
scaledata <- t(scale(t(z))) # Centers and scales data.
hc <- hclust(as.dist(1-cor(scaledata, method="spearman")), method="complete") # Clusters columns by Spearman correlation.
 TreeC = as.dendrogram(hc, method="average")
 plot(TreeC,
 main = "Sample Clustering",
 ylab = "Height")

In fact the cpm function will use the size factors (TMM) that were calculated with calcNormFactors. In the absence of the size factors cpm would do a naive per-million scaling so only correct for differences in library size. The power of the size factors is to also correct for library composition. I recommend this video for details. In short: What you do is fine. I would suggest though that you use the log = TRUE option for cpm as one typically performs downstream applications (even when using the Z-score) on the log scale. Or alternatively t(scale(t(log2(z+1)))).

Thank you for your comment. Could you upload the video link?? I cannot see...

But TMM also conduct the library size normalization, right??
Can we say this is a kind of redundancy?? What is a role of cpm transformation after TMM??
https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html

No it is not redundant. The normalization in edgeR works by first scaling for equal library size and then additionally applying the size factors to correct for library composition. The calcNormFactors calculate the factors, the cpm function applies it to the per-million scaled counts.
I refreshed the link: https://www.youtube.com/watch?v=Wdt6jdi-NQo


text:
WGCNA co-expression network analysis with less than 20 samples
I have two cancer subtypes data. Subtype A is 14 samples and Subtype B is 23 samples. I'm interested in identifying the functions of some LncRNAs in the Subtype A group. For this I'm using all protein-coding genes and interested LncRNAs expression data (normalised counts).
From the WGCNA tutorial, I see that minimum number of samples needs to be 20. So, I guess using only Subtype A 14 samples for the analysis is not possible.
Can I use both Subtype A 14 samples and Subtype B 23 samples together for the analysis? In case if I use both, once I get the modules how do I know that the genes in modules are related to Subtype A or Subtype B? The pathways I get from the interested module is related to Subtype a or Subtype B? How this can be identified?

It's hard to say without knowing how different the subtypes are. If you do a common WGCNA, you may find modules related to the differences between A and B as well as to drivers of expression variation in A and B (common and separate). Since A has fewer samples, the variation will have to be stronger to be reflected by WGCNA modules.
You can run the WGCNA, find modules that vary the most in A samples (these may or may not also vary in B samples), then see whether LncRNAs belong to these modules and have a high kME (module membership).

thanks a lot Peter


text:
which is the right prosite pattern?
I would like to ask if the right Prosite pattern for this multiple alignment:

Is it this:
A-T-[AT]-G-x-C-[AGC]-C-x(1,4)-A

or this:
A-T-[AT]-G-x-C-[AGC]-C-x(1,3)-[ATGC]-A

Actually, I would like to ask if x always implies a gap, and if it is proper to use it when any symbol but not gap can be in the sequence. Thank you!

If nucleotides are acceptable in prosite format because I thought it was protein motifs, the most accurate would be,
A-T-[AT]-G-x-C-[AGC]-C-[GCA](0,1,3)-x-A

The key point is that,

[GCA](0,1,3) of x(0, 1, 3) means zero/ nucleotides (3 indels) or one/ G or C or A nucleotide, and 2 indels or three/ G or C or A nucleotides.
x is any nucleotide i.e. A or G or T or C - it doesn't mean an indel. x(0, 1) is the best representation of a single indel (you call them gaps), here meaning any nucleotide may or may not be present at a given site.

Example from Prosite documentation:

<A-x-[ST](2)-x(0,1)-V
This pattern, which must be in the N-terminal
of the sequence ('<'), is translated as:
Ala-any-[Ser or Thr]-[Ser or Thr]-(any or none)-Val


text:
How to tell if library generation for RNA-seq experiments are stranded vs not stranded
I would like to use RNA-seq data from the NIH GEO database. How can I tell if library generation for the RNA-seq experiments are stranded vs not stranded?

If you already have a reference transcriptome assembly, then you may run the Salmon tool for quantification using your RNA-seq samples. It has an automated library type inference function (https://salmon.readthedocs.io/en/latest/salmon.html#what-s-this-libtype). It might help you to determine whether the libs are stranded or not.

Unfortunately, the stranded nature of cDNA sequencing is not typically included in the metadata attached to a project, which means that downloading a small subset of the data and checking for strand consistency (as suggested by the other answer) is often the quickest way to work that out (when that information isn't included in the methods section of the associated paper). This knowledge is not particularly useful for casually browsing projects to find a stranded protocol that fits a particular study specification.
Sometimes you get lucky, but not often. GEO RNASeq datasets usually have links to the SRA accessions, and those accessions usually include annotation about the library type. For example, consider this dataset on GEO:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE151535
It links to SRP265460 on SRA. Following the link for the first run for this project, I get to here:
https://www.ncbi.nlm.nih.gov/sra/SRX8436504[accn]
This includes the following metadata about the run:

Instrument: NextSeq 500
Strategy: RNA-Seq
Source: TRANSCRIPTOMIC
Selection: cDNA
Layout: SINGLE
Construction protocol: Replicates of oropharyngeal explants at different stages of development (37/38, 39, 41) were immersed in RNALater. Samples were homogenized with a Dounce homogenizer in RLT buffer containing ß-mercaptoethanol (Qiagen). The homogenate was then processed using an RNAEasy kit (Qiagen, Hilden, Germany). Replicate samples for each stage were analyzed for RNA quantity and quality on a Nanodrop 8000 spectrophotometer. RNA libraries were prepared for sequencing using standard Illumina protocols

But, as you can see, there's nothing about whether its stranded or not: Illumina has both stranded and unstranded "standard" protocols.


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:
Residue coevoution from position wise residue preferences rather than multiple sequence alignments
Hello experts in residue coevolution analysis of proteins,
I am looking for a way to calculate residue coevolution (coupling scores) between all pairs of residues of a protein.
The issue is that rather than having a Multiple Sequence Alignment (MSA) which is used as an input by the relevant tools (e.g. evcouplings), I have position wise preferences for 20 amino acids for the protein.
For the protein with length L, I have a 20*L matrix of position-wise preferences for amino acids.
I would like to use this matrix as a direct input for coevolution analysis.
Would it be possible to calculate the residue coevolution direcly from a residue preference matrix?
Any recommendation for a method of choice or a software tool would be helpful.

Coevolution requires you to know when you have residue X at position n you will preferentially have Y at m. Namely the covariance of n with m. If you have a PSSM matrix like you have, you cannot infer the covariance. You would need a L*L*20*20 tensor.

However, your PSSM-like matrix contains useful information, such as what is core and what is surface. But I am not sure what tool uses it. Have a gander at https://predictioncenter.org/casp14/index.cgi


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

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

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


text:
Step-by-Step Construction of Gene Co-expression Networks from High-Throughput Arabidopsis RNA Sequencing Data
I am following the tutorial here to learn how to process raw RNA-seq data and get gene read counts. I've had no issues until I got to step 3.9, "Assignment of RNA-Seq Reads to Genes." Here, the authors run the command
fc0 <- featureCounts(sam.list, annot.ext =
"Ara-port11_GFF3_genes_transposons.201606.gtf", allowMultiOverlap = T, isPairedEnd = F,
nthreads = 2, strandSpecific = 0, isGTFAnnotationFile = T)

I ran the same command except I set isGTFAnnotationFile = F because I couldn't find the gtf file, so I used the gff file TAIR10_GFF3_genes_transposons.gff instead, downloaded from here (they mentioned that using gtf only speeds up the process but that you can still do this using gff). I'm not sure if this is what is causing the issue, but the error I get is:
Error: Line 1 contains a format error. The expected annotation format is SAF.
Error in featureCounts(sam.list, annot.ext = "TAIR10_GFF3_genes_transposons.gff",  :
  No counts were generated.

Thanks for answering the question @AnIgnorantWanderer please consider "Accepting" your answer below because it helps with site statistics.

Answer from @an-ignorant-wanderer, converted from comment:
I ended up using the annotation file from https://support.illumina.com/sequencing/sequencing_software/igenome.html and it worked.
[Please update/improve this answer if possible]