text:
Is removing samples based on clustering for downstream analysis a right choice?
I'm using TCGA Lung cancer data. I'm interested in doing differential analysis between Lung vs Normal. Before DEA, to check the distance between each pair of samples I plotted an MDS plot:

In this, I see some Tumor samples are clustered with Normals. As this is a two dimensional plot, it is a little hard to say.
So, I made a clustering heatmap with the top 10% genes:
 (Where brown color is Normal and black color is Tumor). I also see some Normals clustering with Tumor here, and some Tumor are with the Normals.
My question is: Should I exclude those tumor samples which are clustering with the normals for differential analysis? Is it a good idea to remove samples based on clustering?
Note: I didn't see any technical bias of those samples

I would be very hesitant to blindly exclude those samples based on the clustering. Check to see if the clusters actually denote some sort of batch effect, since it's not like all of the TCGA datasets were processed at the same place or time. Check for clinical covariates too. Also, do a PCA and check the genes with high loading on the relevant PC. If it appears that the genes are indicators for cancerous state, then it's likely the samples are just mislabeled and you can reasonably exclude them. Otherwise keep them in so your results better represent the underlying biological variability.
BTW, your heatmap and MDS plots are what I'd expect from a cancer dataset. You occasionally see a tumour with the normal samples and vice versa, since the problematic genes aren't always in the top 10% and there is likely some heterogeneity in the cancer samples.

How do I check the batch effect of the samples? can I check that with edgeR or Deseq2?

Go through the meta information and see if there are other things (sequencing data, processing center, etc.) that correlate with the clusters.


text:
How to find correlation between two specific genes in same dataset?
I would like to plot the correlation between two specific genes in my data. I have a matrix with genes in rows and samples in columns, with read counts data. If I want to look at the correlation between 2 genes
1) do I need to convert the counts to some other units?
2) Is normalization needed?
I currently have counts returned by featureCounts:
          Sample1  Sample2  Sample3  Sample4  Sample5

ESAM         803    338       528       841      712
ESCO1        1594   1401      847       392      611
ESCO2         274   778       288       77       204
ESD          6192   3502      2013      1573     1808
ESF1         1356   1497      1011      372      701
ESM1         398    58        114       582      135

I tried this way:
I transposed the above data making genes as columns and sample names as rows. I removed the rows. Lets say the genes as columns is in dataframe "df".
df.mat <- cor(df)
colnames(df.mat) <- NULL
rownames(df.mat) <- NULL
heatmap.2(df.mat, trace="none")

I have an error:
Error in .External.graphics(C_layout, num.rows, num.cols, mat, as.integer(num.figures),  :
  invalid graphics state

I want to look at the correlation between ESCO1 and ESCO2. Can anyone tell me how to do that and plot with a small example.

I suppose in an idealized world you'd use something like salmon to get TPMs, which can account for differences in effective gene length by sample. However that's quite likely to be overkill.
Not necessary, you're using a with-sample ratio, which wouldn't be changed by normalization.

Here is a small example in R, assuming your matrix or data frame is named m:
first = which(row.names(m) == "ESCO1")
second = which(row.names(m) == "ESCO2")
plot(m[first,], m[second,])

or cor(m[first, ], m[second,]). You'll want a rank-based correlation.

When I did "cor(m[first, ], m[second,])" I see column names as both rows and columns and showing "NA" for all.

Stop removing the column and row names.

don't know what's wrong. the dataframe "df" is with column names "gene", sample1, sample2, sample3 etc...I made the first column "gene" into rows. then I got the "ESCO1" into first and other into second. Then I did cor(df[first, ], df[second,]) and didn't work.

"It didn't work" isn't ever a useful thing to write unless immediately followed by either the error message or problematic output you received. I presume your dataframe still has the gene column, which would need to be removed, or you need to put things through c() to get a proper vector.

@beginner when asking this sort of thing, you need to give us an example that reproduces your issue which we can run on our own machines. You are running some command with a wrong parameter, but since you don't show us the code you use, we can't help you find the mistake. Please see https://stackoverflow.com/help/mcve

@terdon sorry I misunderstood. but I tried in an other way and answered my own question. Please check that and let me know. thank you

@DevonRyan I tried in the following way. Please check that. answered my own question. thank you

You can represent, and calculate, the correlations using psych:
library(psych)
pairs.panels(df,df,
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE,xlim=c(0,500),ylim=c(0,500))

I'm answering my own question. I tried this way.
df <- read.csv("df.csv")

Samples ESAM    ESCO1   ESCO2   ESD ESF1    ESM1
Sample1 803     1594    274    6192 1356    398
Sample2 338     1401    778    3502 1497    58
Sample3 528      847    288    2013 1011    114
Sample4 841      392    77     1573  372    582
Sample5 712      611    204    1808  701    135

df2 <- data.frame(df[,-1],row.names = df[,1])

cormat <- round(cor(df2,use="pairwise.complete.obs"),2)
head(cormat)

library(ggcorrplot)
pdf("eg.pdf")
ggcorrplot(cormat,hc.order = TRUE,
           outline.col = "white",
           ggtheme = ggplot2::theme_gray,
           colors = c("#6D9EC1", "white", "#E46726"))
dev.off()

This does not answer your own question, you ask for a plot with the correlation between 2 genes, but you show a correlation plot between all your genes...

yes, this shows correlation between all the genes but also shows the genes I'm interested in. I don't know how to plot for a specific genes so , I did this.

Down vote from me, you ask for something else in your question (plot between two genes). Impossible to 'guess' what you are really meaning. Also from your answer it is not clear if you have to normalize or not.

source: Stack Exchange

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

sort: votes


text:
Ensembl id to GeneSymbol with biomart
I have 3224 Ensembl id's as rownames in a dataframe "G". To convert Ensembl ids into Genesymbols I used biomart like following.
library('biomaRt')
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- rownames(G)
G <-G[,-6]
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=genes,mart= mart)

Now in G_list I can see only 3200 ensembl ids showing Genesymbols / No Gene_symbols. Why the other 24 ensembl ids are not seen in G_list? If there are no gene_symbol for those 24 ensembl ids it should atleast show "-"
Examples of problematic IDs are: ENSG00000257061, ENSG00000255778, ENSG00000267268. These are not at all shown in G_list (biomaRt). So, I gave them in biodbnet, which seems to handle them.
what is the problem here?

I used gProfileR where there are no problematic IDs with that. It converted ENSEMBL IDs to Gene_symbols and made my work easier with GO analysis too. There is an R package and also API.

could you post a code snippet? or is it trivial?

Here is the reference manual for gProfileR [https://cran.r-project.org/web/packages/gProfileR/gProfileR.pdf]

It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have - as symbols don't have associated symbols, but are present in the database.
To use an archived version in biomart:
mart = useDataset("hsapiens_gene_ensembl", useEnsembl(biomart="ensembl", version=84))

That's an example for release 84.

biomaRt will return NA for genes that are in the database, but don't have symbols.

@IanSudbery Yeah, I'm assuming that OP replaced NA with - at some point.

As you said, I thought the same that the problematic ID's are existed in GRCh37. And tried with that too. And among 3224, I found only 3089 in the results.

Those IDs existed in the GRCh38 release 78 annotation (just checked a local GTF), so perhaps they're still GRCh38, but just an older release of it.

I actually did it in this way. install.packages("devtools")
devtools::install_github("stephenturner/annotables")
library(dplyr)
library(annotables)
grch38

So ENSG00000257061 was last seen in ENSEMBL 84 (you can look this up by searching on the ENSEMBL website). You need to search against the exact same version as you got your gene list from. Replace useMart with useEnsembl(version=X) where X is the version of ensembl you used to generate the gene list.

So...can anyone of you please tell how can I use the latest version in the above code? or Is there any R package for biodbnet?

@user3351523 I've updated my answer. The biodbnet question should be a separate post.

Just by-the-by, I'm going to guess that any ID that has been retired is unlikely to map to a HGNC symbol.

@DevonRyan With the above mentioned version 84 I got 3220 now. Still 4 ids missing.

Try version 83 then.

@user3351523, you have to use the same ensembl version as was used for generating your genes object, or you'll never get 100% conversion.

I used "gProfileR". There are no problematic IDs with that. It converted ENSEMBL IDs to Gene_symbols and made my work easier with GO analysis too. There is an R package and also API. Thank you all for the help !!

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


text:
How to identify gene expression signatures from gene expression data?
I have TCGA gene expression data. I'm interested in identifying gene expression signatures using the data.
I would like to know whether there are any tools or R packages for identifying gene signatures.
How are gene signatures different from GSEA?

@gringer, usually these is the wording used to denote genes that their expression is characteristic of a specific group of cells (cell lines or not), I read it also as marker genes (not in the sense that are used because one can mark them in the surface or in a specific region of the cell but because the expression pattern is specific of these cells). But user335123 Could you confirm that we are talking about the same meaning?

@Llopis yes what said was right !!

Why do you want to do this? I have never seen any reason to generate new gene signatures/gene sets from TCGA, since they already exist for humans. Or are you trying to create signatures for specific cancer types?

BTW, if the (more standard) definition of a gene signature from @Llopis is correct then please edit your post accordingly.

@Devon, do you mean that exists other gene signatures for humans (in a database)? Could you please provide me a link or a reference? User I'will answer later

@Llopis The question needs to be appropriately clarified by user3351523 first, since whether a database exists is dependent upon that (and anyway I'm not putting it in a comment that I'm going to delete).

I don't understand: my comment got deleted and now the conformation that what I rephrased is isolated (you where refering to that comment that could be deleted?). Are you looking for clarification on if it either for human or if it is what I had written in my first comment?

@Llopis Oops, restored the comments. They got flagged for deletion and I didn't check thoroughly enough. OP needs to correctly edit into their question what the exact goal is, since it seems that they didn't even define gene signature properly in the original post.

GAH, the incorrect definition wasn't from OP!

@DevonRyan I'm interested in generating new signatures for lung cancer histological subtypes.

It was a comment from OP, which I flagged for deletion after adding it to the question. The question is not clarified enough as it is to give a useful answer.

There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a gene signature:
Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down regulated genes and defines them as the signature of the group of cells. That's how the gene sets of the MSigDB which are not from (or derived from) pathway databases or gene ontologies are defined. There is no gold threshold of a fold change or a specific number of genes, it (highly) depends on the experiment you run and how good are the isolation of the cells.
Tools: DESeq2, limma, edgeR,...
As suggested by gringer, another way is to find which genes are expressed together (independently of how do they relate to other groups). This method relays on having a good enough number of samples from the same condition to be able to identify which genes are expressed the same way in those cells. If these genes are not expressed that way in other cell lines, then you have your signature for this group. Again, there isn't a golden method, different types of co-expression networks can be build and depends on how specific do you want to be or which is the data provided.
Tools: WGCNA, diffcoexp, TReNA,...
This is different from a gene set enrichment analysis1, which consists on evaluating if a group of genes (usually a gene set, a pathway or genes from a specific gene ontology) are not randomly sorted in a sorted list of genes (usually this sorted list is the fold-change of a differential analysis).
There are several methods developed depending on the null hypothesis, the type of sub-method and what do they take into account...
Tools: GSEA, PLAGUE, GSVA, limma, fgsea,...
TLDR: To perform a gene set enrichment you don't need a gene signature but a gene set (that can be a gene signature or a pathway or a random group).
1: Note that gene set enrichment analyisis is different from gene enrichment analysis. For a good review on the methods I recommend this article

If you're interested in an unbiased approach you could have a look at Gene Correlation Network Analysis. WGCNA is a fairly old R package, with a somewhat confusing workflow and results that can be difficult to interpret. It also needs a large sample set, preferably from fairly different biological conditions, to produce meaningful results.
Asking how "gene expression signatures" differ from GSEA is like asking how surface textures differ from trademarks. Gene set enrichment analysis is for looking at gene rankings and seeing how they compare to a known / existing gene set. A gene set is considered "enriched" if the majority of the gene set appears in high-ranked genes (or alternatively in low-ranked genes), rather than a more even distribution of the set across the ranked gene list. Most commonly people will rank genes based on expression (or fold change), but the method is agnostic to the nature of the ranking.

WGCNA does need a large sample set but from the same biological condition in order to find the meaningful correlations between genes, otherwise the correlations are not good enough to infer that there is a relationship between genes

What does large sample set mean here? how many samples you mean?

I'll start the ball rolling with a suggestion of 40 samples: 40 cancer subtypes that are a good representation of the biological variation that you'll be testing on

Gene Signatures are sets of genes to report a consensus signal or function. They are not defined in a particular way but can be defined by any chosen group depending on the biological function that you are interested in. Typically these are defined by differential gene expression analysis or principal component analysis, although other approaches such as network analysis or time-course data could also be used.
There are plenty of examples of sets of genes published as gene signatures (especially in the field of cancer research). Those that have been cross-validated in other datasets or replicated by other research groups are ideal. If you wish to generate such a “signature” for a particular biological function, you could follow the approaches demonstrated in these published signatures (especially those that have been widely used and cited).
A gene expression signature is used to generate a consensus expression signal for these groups of genes (representing the given biological function). A widely accepted way to do this across samples of bulk RNASeq (or microarray data) is a “meta-gene” from the first eigenvector (v[,1]) of the singular value decomposition (SVD). This dimension reduction identifies a major contributor to variation across samples (in these genes) without negatively correlated genes averaging out. Dimension reduction can also be achieved by PCA and tSNE, depending on the requirements of your analysis.

Hi- I'm interested in the “meta-gene” approach using SVD. Do you have any recommended reference at hand describing it? With the occasion, another approach to identify "meta-genes" is based on non-nagative matrix factorization.

It’s described rather well here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1261155/#!po=65.1786

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

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


text:
Mathematical notation for formulating a rank score
This is a cross-post from Mathematics forum. As no one has answered it yet, I fIgured post here as well.
I would like to describe a transform I used to rank my data points. I have recoded my variable columns from $10$ measurements with different range to $n$ (for each measurement) where $n \in \{0,0.5,1\}$.
For columns where measurements were expressed as categorical values I simply recoded the categorical value to either $0$, $0.5$ or $1$ and for continuous variables I recoded the lower quartile as $0$, interquartile as $0.5$ and upper quartile as $1$ and finally I summed the recoded values to produce a single score $\sum n_{i1,..i10} $ for each row.
I am trying to write my method and I'd like to have it in mathematical notation. I was wondering if anyone could help me with this!

Maybe try https://stats.stackexchange.com/

What do you mean by "have it in mathematical notation"? Can you give an example of something you consider "in mathematical notation"?

Also, it would be good to include the link to the cross-post

Hi bli, Thank you for your comment. By mathematical notation I mean mathematical representation of my transform. I am including the link as you suggested!

Your question is not very clear to me. So can you please try to explain it with some example. Then we may able to help you.

Have you looked at the dataframe I uploaded in Mathematics forum?

Yes, I did. But still it not clear.

I think the problem with your question is two-fold. 1) you don't provide any reproducible data so that we can grasp the details of how you calculate this "rank" score. See this post for an example how you would make a reproducible example in R to get an idea of what we mean by "reproducible example" 2) there is not likely to be a mathematical notation for what you have done, because there is no mathematical formula to convert strings ("I", "II", "III") into numbers. One way that you could potentially solv

I think that "mathematical notation" does not necessary have to consist in "formulas". One could probably define a bijection between a set of strings and a set of numbers, in proper mathematical notation, and that would validly describe a conversion procedure.

Yeah you are probably right that there is some valid way of describing the rank score. I imagine in this case it would be something more akin to a key than notation/formula, but if you have an idea how to describe it why don't you add an answer?

As you pointed out, it is difficult to grasp the details of what the OP wants. Besides, it's been a long time since I last wrote maths formally.

Here is my solution to the problem. I am posting it here in case anyone else come up with the same idea but did not know how to formulate it in mathematical notation!

$$ \psi_i = \sum_{i=1}^9 S_i $$
where $S_i$ is the score function and defined as:
$$
S_i = \cases{
  0 & \text{if $\theta(x_i) < Q_1$ or $\theta(x_i) = $ benign/neutral} \\
  0.5 & \text{if $Q_1 \leq \theta(x_i) < Q_3$ or $\theta(x_i) = $ possibly damaging/uncertain} \\
  1 & \text{if $\theta(x_i) \geq Q_3$ or $\theta(x_i) = $ damaging}
}
$$
The $\theta(x_i)$ is the pathogenicity or conservation score for variant $x$ as defined by model $i$ and $Q$ denotes the quartile range for scores from M-CAP, CADD, GERP and Phylop models.

If Q are quartiles, then the last line for damaging shouldn't be $>Q_3$?

@Konrad and Reza : Can you look into this problem and provide some suggestions. https://math.stackexchange.com/questions/2716192/how-to-convert-a-detailed-mathematical-statistical-compuational-process-into-ma


text:
Why total RNA-seq usually yields low mapping rate?
Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)?
Doesn't the genome fasta file contain ribosomal RNAs (which are abundantly expected in the total RNA-seq libraries) as well?
I don't have the numbers for human right now, but I recall previously having low mapping rates for human data as well. What triggered me at the moment is actually when I mapped zebrafish data. I got around 60% mapping to the genome and 45% to the transcriptome with one dataset and 36% - 35% with another dataset (note that these are from 2 different studies and both are total RNA-seq, mapping was done using STAR and Salmon, which have their own methods).

It can also happen if you map to a transcriptome sequence instead of the genome sequence because there can be a lot of intron sequence from unprocessed mRNA and from unannotated non-coding RNAs

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (this behavior can be changed with the --outFilterMultimapNmax option).
To confirm if this is actually the case, you could check the number of multi-mapping reads in the aligner's log files.

I tried to increase the number of allowed multimapping reads and indeed the number of reads mapped to multiple loci has increased, albeit slightly (6.92% to 7.50%). Most of the non-mapped are categorised into "too short" (in STAR). Do you have any ideas what are these (well, it refers back to my initial question, since this happen in Salmon as well)?

Are you trimming adapters from the reads prior to mapping? That might explain the many "too short" reads. Also, have a look at this issue on STAR's Github repo: https://github.com/alexdobin/STAR/issues/169

@kaka01 If accounting for multi-mapping doesn’t solve your problem then there may simply be something wrong with your data: on high quality data sets, mapping total RNA to a genomic reference should typically yield >80% mapped reads. Many “too short” mappings may be indicative of RNA fragments. Did you perform size selection before sequencing?

Tom: No trimming was done prior to mapping. @Konrad The libraries weren't size selected. Your explanation makes sense, aligning the unmapped reads to the collection of rRNAs yielded high mapping rate (which is a sign that those are really rRNAs). Still not sure why STAR classified them as too short, maybe I just have to play around more with the parameters.

@kaka01 “Too short” here simply means that STAR is unable to align the read with high quality: either because the initial read (after trimming) is so short that it could match the reference virtually anywhere (< ~14 bases for something like the human genome) so we have low confidence in the correct origin. Or because STAR, when running with --alignEndsType Local (which is the default) is only able to match a small part of the read (which then leads back to the first problem).

@kaka01 Since you didn't do size selection, it's possible that you have a high number of short fragments, and you end up sequencing the adapter. You could try to trim the adapter and map the reads again. This might help to mitigate the issue because in STAR the minimum mapped length if a fraction of the total read length.

@KonradRudolph in you answer you say that mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. I wasn't aware of this issue, do you know why this is the case? (sorry for commenting here, I don't have enough reputation to comment you answer directly)

@Tom Let’s talk offline (are you at the RNA club tonight? I’ll attend.) In short, I have no confirmation but it’s what we observe, and I don’t find it implausible given the chemistry.

In particular for mouse, this can happen if the ribo-depletion wasn't terribly efficient, since there's no Rn45s sequence in the reference genome. That, combined with numerous copies of tRNAs, 5S rRNAs, etc. causing issues with multimappers (see the answer from @Tom) can heavily decrease alignment rates.
Note that the human reference genome contains a few copies of the 45S (e.g., on GL000220.1 and chr21), so these reads will all multimap (possibly too many times). I don't know about the zebrafish reference genome.

When sequencing RNA from freshly harvested tissues under proper conditions, you should generally expect > 50% mapped reads. In fact, everything < 80% would usually raise concerns.
From your description (in the question and comments) it sounds as if your samples are potentially degraded and therefore saturated with short RNA fragments, either because you’ve got old tissue, you’re looking at environmental samples, or from improper handling.
And since you’re not performing any enrichment, the sequencing data will then be highly saturated with degraded RNA, too. The reason for this is that RNA-seq is a (hopefully somewhat uniform) sampling of your sample aliquot: everything in it is proportionally represented in the RNA-seq data. That’s why enrichment steps are important, to increase the proportion of whatever fraction we’re interested in.
In particular,

Total RNA is mostly composed of genes with many copies: mostly rRNA and (much less) tRNA. Hence the importance of ribosome depletion or poly(A) selection.
RNA in general is unstable. Depending on the sample collection (non-tissue sample origin? age of the tissue before harvesting?), a large fraction of the RNA may have degraded into minuscule fragments, which are further shortened by the library preparation. Size selection would get rid of these undesired fragments.

In my experience, the second issue is exacerbated on certain sequencers (I’ve seen it on an Illumina HiSeq 1500) because a mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. So even if you have a (small but still present) fraction of long RNA fragments, the sequencer may not be able to synthesise their whole length efficiently. The effect is that there are proportionally more short reads than there were short fragments in the sample.
You can verify this by calculating the insert size distribution of your reads after adapter trimming (e.g. using picard, or by simply tallying read lengths1). Small reads (< ~14 nt) from degradation fragments are essentially unmappable, since they are too short to have a specific identity. Even if we allowed the mapper to align them their coordinates would be essentially random.

1 This works:
awk 'NR % 4 == 2 {c[length($0)]++} END {for (i in c) print i, c[i]}' in.fastq

You can't assume that the rRNA genes are in the genome. There many copies of them and not all are placed on chromosomes. We've had this issue with human and mouse data for total RNA preps. You need to ensure you map against the whole genome not just the chromosomes.

I 'm mapping to the whole genome indeed. But anyways, assuming that you are mapping against the whole genome, did you get significantly higher mapping rate?


text:
Counting number of possible alignments between two DNA sequences using python
I was looking BioPython to calculate the number of possible alignments between two sequences. Let's say there are two input sequences with the length of m and n. The program will count the possible number of alignments based on the length of two sequences.
Is it possible to do that with BioPython? If not, how can I do that?

Which kind of sequences DNA, RNA or proteins? How have you tried to do so? What do you expect to learn from it?

And the third question? Why do you need to know? (By the way, if you don't try why do you expect random people on the internet to help you?, Not that someone can't but it is harder to provide help if one doesn't see research effort)

Seems like a homework assignment where they are asking for the theoretical number of alignments.  even down to the typical n and m nomenclature.

You should define what you mean by "possible alignment". Biopython will only help you use tools that chose among the possible alignments, based on certain criteria.

A closed-form solution is offered in An exact formula for the number of alignments
between two DNA sequences by Torres, Cabada, and Nieto:
$$f(m,n)=\sum_{k=0}^{min(m,n)}2^{k}\binom{m}{k}\binom{n}{k}$$
If this solution seems reasonable, you could calculate this without BioPython, but with a simple xrange loop, the power operator, and scipy.special.binom:
#!/usr/bin/env python

import sys
import scipy.special

def alignments(m, n):
    if m < 0 or n < 0:
        raise ValueError('m and n should be non-negative')
    s = 0
    b = min(m, n) + 1
    for k in xrange(0, b):
        s += (2 ** k) * scipy.special.binom(m, k) * scipy.special.binom(n, k)
    return s

def main():
    m = 4
    n = 2
    a = alignments(m, n)
    assert (m == 4 and n == 2 and a == 41), "something went wrong!"
    sys.stdout.write("alignments(%d, %d) = %d\n" % (m, n, a))

if __name__ == "__main__":
    main()

I don't know what your m and n are, but keep in mind that this solution grows quickly. Also, I did some digging and scipy approximates the calculation of the binomial coefficient, which is probably going to be an issue for large m and n.
If you're working with large values, you may want to look into using the Python long type and consider the accuracy of the answer. The Stack Overflow link offers alternatives to scipy that generate an exact calculation of the binomial coefficient, at the cost of a longer runtime.


text:
How do I export my biom object to a file in biom format?
I'm using the python API for the biom-format package. I read in two files in biom format and merged them, like so:
In [1]: import biom

In [2]: a = biom.load_table('./table-a.biom')

In [3]: b = biom.load_table('./table-b.biom')

In [4]: a
Out[4]: 51116 x 46 <class 'biom.table.Table'> with 185077 nonzero entries (7% dense)

In [5]: b
Out[5]: 37479 x 26 <class 'biom.table.Table'> with 100124 nonzero entries (10% dense)

In [6]: c = a.merge(b)

In [7]: c
Out[7]: 67022 x 47 <class 'biom.table.Table'> with 254555 nonzero entries (8% dense)

How do I now export this biom object (c) as a biom file? There are a few methods available, eg, .to_tsv(), but I don't see a way to save this object in biom format.

Are you using the biomformat package of bioconductor through Python? Or the biom-format package? Anyway in R the relevant function is write_biom so it might be the same in the python package

You can use the functions to_hdf5 and to_json from the biom-format package:
to_hdf5
with biom.util.biom_open('output.biom', 'w') as f:
    merged_table.to_hdf5(f, "example")

to_json
with open('output.txt', 'w') as f:
    f.write(merged_table.to_json("example"))

There is a function in PiCRUST you can use
from picrust.util import write_biom_table

write_biom_table(c,'table-c.biom')

Note, however, that depending if you want json or hdf5 format and if your metadata are complex, you may need to use additional arguments. These are the arguments of the function:
def write_biom_table(biom_table, biom_table_fp, compress=True,
                 write_hdf5=HAVE_H5PY, format_fs=None):
"""Writes a BIOM table to the specified filepath

Parameters
----------
biom_table : biom.Table
    The table object to write out
biom_table_fp : str
    The path to the output file
compress : bool, optional
    Defaults to ``True``. If True, built-in compression on the output HDF5
    file will be enabled. This option is only relevant if ``write_hdf5`` is
    ``True``.
write_hdf5 : bool, optional
    Defaults to ``True`` if H5PY is installed and to ``False`` if H5PY is
    not installed. If ``True`` the output biom table will be written as an
    HDF5 binary file, otherwise it will be a JSON string.
format_fs : dict, optional
    Formatting functions to be passed to `Table.to_hdf5`

Notes
-----
This code was adapted from QIIME 1.9
"""

I have some problems to perform a simple operation and write the output of biom data generated with picrust, it might be helpful as an example of how to use the arguments.


text:
How can I calculate the silent mutation for each position in codon?
I am trying to find the number of silent mutations in for each position in codon.
In the homework, I am asked:

For the codons and for each position, count the number of silent mutation. There are 64 codons, for each position there are 3 mutations leading to a change at the DNA level. There are thus 192 possible mutations for each of the three positions of the codons. For each position, print the number of silent mutations.

I understand stand silent mutation is the mutation that doesn't lead to amino acid change. For example, GCC can become GCT but it is still Thr
I am not sure what does "calculating silent mutation for each position" mean? I mean it is easy to find silent mutation for each codon like the example I gave above, but couldn't really understand how to calculate the silent mutation for each position.
Should I change letters (nucleotides) one by one and see if it causing a silent mutation?

Calculate in this context would just mean determine. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.
BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. Terming things like this as "silent polymorphism" or "silent change" or something along those lines is preferable.

but he wants me to find number of silent mutations for each position of codon. Lets say codon is GCT and wants me to count the number of silent mutation for the codons and for each position.

You can efficiently get that from the method I described.

@snnlaynnkrdsm try and implement this. Once you have something written down, if you're stuck with a specific aspect of your program, you can ask again here. But first try it out! We won't just write the program for you, but we're happy to help with specific problems you face.

Curious why you think 'mutation' is a loaded term? To me, 'silent' is the loaded part, since synonymous changes may have fitness effects due to codon usage bias or regulatory element overlapping coding sequences

Patients (and scientists in my experience) tend to over interpret what a "mutation" is. Things like synonymous polymorphisms were less likely to be conveyed to non-professionals, but yes the "silent" terminology is also misleading, though it's not going to make some random patient think they're going to be the next Magneto.


text:
Which secondary structure prediction algorithms don't query the PDB/DSSP for known homologous sequences?
There are several secondary structure prediction algorithms. I am interested in the algorithms which predict secondary structure without accessing the PDB to cheat by using the DSSP or STRIDE entry for known homologous structures.
Since DSSP and STRIDE are not performing predictions, rather they analyse the PDB structure to determine whether each residue is a helix, strand or coil.  This in turn means that, the 'so called' secondary structure prediction algorithm, also, is not performing predictions (except, of course, in the cases where no homologous structure/sequence is found in the PDB).
Which secondary structure prediction algorithms don't cheat by using DSSP or STRIDE for homologous structures?
Here is a list of secondary structure prediction algorithms, which are accessible via a consensus secondary structure prediction server called @NPS (currently down: https://npsa-prabi.ibcp.fr/NPSA/npsa_seccons.html)
DPM
DSC
GOR1
GOR3
HNN
MPSA
PHD
PREDA
SOPM

Have you tried googling for ab initio secondary structure prediction?

@DevonRyan It is a good idea, and doing so does provide some useful information.  However, ab initio modelling is not what I am interested in.  Use of knowledge is fine, only it should be generalised, not a specific query to the PDB and DSSP for a specific sequence - in such case, no prediction is performed, only cheating i.e. an ~100% accurate answer.

If using knowledge is fine the querying PDB and such when possible is in no way cheating.

Yes it is cheating.  Not performing a prediction = cheating.

You have a warped definition of cheating if it includes using the most reliable information. Since you implicitly don't want the most likely results, what is your end goal with all of this?

A "predictor" not making a "prediction" is cheating.  I am conducting some research into secondary structure predictions comparing their predictive results.


text:
Intron retention events for a condition sample with 3 replicates
I have 3 biological replicates of RNA seq data for a particular condition. I want to find out intron retention events from those biological replicates for a given condition. There is no comparison I am doing here. Just from the three replicates, I am trying to find out what introns are retained in my sample. There is only one condition here. Tools like rMATS require 2 conditions. Is there any tool I can use to find out intron retention events given that I only have one condition (with replicates).
Ideally I would like to use a method similar to that used in this paper, though I can't fully grasp their methodology. This would allow me to get p-values for enrichment of specific introns.

All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.
What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.

I second this, but would suggest using SUPPA which handles replicates out of the box and is actively developed. As for deciding which introns are retained, I expect that most biological conditions most of them will have a low PSI (which a nice, easily interpretable metric), so I have in the past plotted the distribution and used some upper bound of the distribution to be on the safe side. That said, some condition to use as a reference would be better.

The method you linked to is interesting, though I suspect that Ian's answer is both easier and generally sufficient to get what you want.
Having said that, the general idea behind the method you linked to is to use DESeq2 to test for differentially expressed introns. To generate the control samples, they're effectively taking the number of reads across all introns in a gene and distributing them over the introns. They weight the distribution such that longer introns receive more reads (they use an effective length for this, but I'm skeptical that this makes that much of a difference). DESeq2 is then giving the p-values you're after. Unless you can get the authors to send you there code, you'll have to write a bit of code to do all of this...but it shouldn't be too bad.
BTW, you can use the GenomicFeatures package from bioconductor, specifically the exonsBy(..., by="gene) follows by gaps().

Hi Devon,
Thanks for your help.The paper compares insilco replicate counts with actual counts and that is what is input for DESeq2. From your explanation, for a gene, they took all the reads which map introns and then lets say we get 1000 reads for all the possible introns for this gene and then for each intron, they give a weightage like 1000/(len of intron). I am looking at GenomicFeatures and I have already created the intron file coordinates which you mentioned at seqanswers. If possible and if you have time can you guide me in writing the above code. Thanks for your help always
Regards

It's not 1000/(len of intron), but rather 1000/(len of intron / sum of all intron lengths in gene).

I only have time to provide vague guidance on the code. You can use findOverlaps() (from GenomicRanges in bioconductor) for the counts. Aside from that it's mostly split() and mendoapply() that you'll need.

so that would give us insilco replicate1 count for a particular intron and similarly insilco replicate2 count for 2nd replicate and then we compare these 2 with actual counts we are getting in rep1 and rep2 via DESeq2

Yup, exactly. Once you have the counts its just down to using DESeq2 as normal.

D'oh! I've actaully done this before (that'll teach me not to look at the links in the question). I used feature counts to count the reads in introns, then a python script to calculated the effective length from mappability tracks downloaded from UCSC, and then used an R script to find the "detained introns".  See these gists: https://gist.github.com/IanSudbery/df92ebf4f84d76027394750a6516db1b , https://gist.github.com/IanSudbery/eb36ba5d59b4455f9cd7b2085447eb3d. I offer no guarantees with these, they were pulled from my archives.

It doesn't work too bad. Better on nuclear RNAseq than cytoplasmic, but it is very sensitive. Finds introns with only a small amount of retention.

IRFinder can be used to detect the IR events and can easily reach your goal.
Here is an example. IRFinder can calculate the IR ratio to detect IR event. mimirna.centenary.org.au/irfinder/example1.html. (You just need run IRFinder for one condition and each replicate will output a folder). Two conditions can be used to detect the differential IR events. If you have only one condition, you can select the IR events from the common events of the three replicates or select IR events exist at least two of three replicates based on the main output files( named IRFinder-IR-dir.txt). The IR events can be selected by filtering the IR ratio, Intron coverage etc

Hi Renhua, this answer doesn't seem to really address the issue. Could you please elaborate how IRFinder helps when there is one condition with replicates? How would that be done with IRFinder? Many thanks for your time


text:
Help Adjusting Bedtools Output in R
I have been asked to write a script in R to edit, using data.frame, this .fa file that is the result from extracting gene sequences from a repeat masker coordinate file on the genome I am working with. Each input looks as such:
>rnd-5_family-5445_Unspecified(-)
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_Unspecified(-)
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATA
AAAATCTTTCTACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTG
TAATTTACAAAATGAATTACTTGAAGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCT
ACCGATGCCAAAAATTGACGATATTCTATTATTTAATCTAAATGTTAGTCTTTACATACCCTCCCCTA
ATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGCTAAATCTAATTCAATAGATT
GATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATTGCATATAA
CTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATG
AAAGTGATC
>rnd-5_family-2614_Unspecified(-)
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATT
AGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCT
ACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAA
ATGAATTACTT
>rnd-4_family-798_Unspecified(-)
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAA
ATAAACCTTATACTGCTAT

This is the repeat masker coordinates that the above were derived from:
>KB824701.1 417 478 rnd-5_family-5445_Unspecified . -
>KB824701.1 587 1072 rnd-5_family-2614_Unspecified . -
>KB824701.1 914 1129 rnd-5_family-2614_Unspecified . -
>KB824701.1 1138 1225 rnd-4_family-798_Unspecified . -

and ideally I would like it to be in the format of this, for example using the first one:
>rnd-5_family-5445_KB824701.1_417_478

I was told to use data.frame in R but I am not sure exactly how this will work? Can someone help please

Hi Noob, what have you done? Are you familiar with R or other shell programs? This seems as homework could you confirm it is so? I wouldn't like to post a solution to get you later accused to plagiarism... or worse, to provide an answer and you just used it without learning from it.

It is not homework, it is actually a project I am being to asked to help work on while also learning to skills so I have a lot of trouble trying to work and learn at the same time. I used a bedtools command getfasta to extract the sequences from a genome I was given repeat masker coordinates for already

Just checking that I understood correctly, you want to change from >KB824701.1 417 478 rnd-5_family-5445_Unspecified . - to >rnd-5_family-5445_KB824701.1_417_478? Or from  >rnd-5_family-5445_Unspecified(-) to >rnd-5_family-5445_KB824701.1_417_478. Are you comfortable with a other solutions not involving R ?

In case an answer would still be relevant after a year and a half, below is a solution with the Biostrings package.
library(Biostrings)

# loading data.table for the fread() function.
# fread() is not only fast but is quite good at guessing how
# the file is formatted
library(data.table)

# readDNAStringSet() function creates a "DNAStringSet" object
my_fasta <- readDNAStringSet("fasta.fa")

> my_fasta
  A DNAStringSet instance of length 4
    width seq                                                                names
[1]    61 TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG      rnd-5_family-5445...
[2]   485 AATCTGAGATTAGATTTATTTATCTGTTATCC...GGCCAATCATTCAAAAGTTATGAAAGTGATC rnd-5_family-2614...
[3]   215 TTTCTAGTTTGTATTTGAGGAAAGAAAATCCT...TAGAAATTGTAATTTACAAAATGAATTACTT rnd-5_family-2614...
[4]    87 ATGCTTAATTTGCTAAGCATTTAATAACCTTT...GTAGACATTAAAATAAACCTTATACTGCTAT rnd-4_family-798_...

# the fasta header is stored at the @NAMES slot of the @ranges slot of the
# DNAStringSet object created above
> my_fasta@ranges@NAMES
[1] "rnd-5_family-5445_Unspecified(-)" "rnd-5_family-2614_Unspecified(-)"
[3] "rnd-5_family-2614_Unspecified(-)" "rnd-4_family-798_Unspecified(-)"

# removing the trailing "_Unpecified..."
my_fasta@ranges@NAMES <- sub("_Unspecified.*",
                              "",
                              my_fasta@ranges@NAMES)

my_coord <- fread("repeat_masker.txt")
> my_coord
           V1   V2   V3                            V4 V5 V6
1: KB824701.1  417  478 rnd-5_family-5445_Unspecified  .  -
2: KB824701.1  587 1072 rnd-5_family-2614_Unspecified  .  -
3: KB824701.1  914 1129 rnd-5_family-2614_Unspecified  .  -
4: KB824701.1 1138 1225  rnd-4_family-798_Unspecified  .  -
#removing the ">" sign

my_coord$V1 <- sub(">",
                   "",
                   my_coord$V1)

# appending first 3 columns of the coordinate file to "NAMES" as requested by the OP
my_fasta@ranges@NAMES <- paste(my_fasta@ranges@NAMES,
                               my_coord$V1,
                               my_coord$V2,
                               my_coord$V3,
                               sep = "_")

# writing the modified DNAStringSet object into a fasta file
writeXStringSet(my_fasta, "edited_fasta.fa")

>rnd-5_family-5445_KB824701.1_417_478
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_KB824701.1_587_1072
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTA
CCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTTGA
AGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCTACCGATGCCAAAAATTGACGATATTCTATTATTTAA
TCTAAATGTTAGTCTTTACATACCCTCCCCTAATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGC
TAAATCTAATTCAATAGATTGATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATT
GCATATAACTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATGAAAG
TGATC
>rnd-5_family-2614_KB824701.1_914_1129
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATTAGATTTATTTAT
CTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTACCACTTTTGATTCATTTTTTAAA
TCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTT
>rnd-4_family-798_KB824701.1_1138_1225
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAAATAAACCTTATA
CTGCTAT


text:
10X cellranger count, [error] The chemistry was unable to be automatically determined
While running cellranger count using the following .slurm file:
#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=48:00:00
#SBATCH --mail-type=end
#SBATCH --mail-user=nv4e@virginia.edu
#SBATCH --output=generate_count_%j.out
#SBATCH --error=generate_count_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

# Run driver program
FASTQS="/scratch/nv4e/kipnis/AHN3LLBGX5/outs/fastq_path/HN3LLBGX5"
TR="/scratch/nv4e/kipnis/refdata-cellranger-mm10-1.2.0"

cellranger count --id=17R_1 --transcriptome=${TR} --fastqs=${FASTQS} --sample=17R --expect-cells=2000

I got the error:
>2018-04-28 20:04:37 [runtime] (progress)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMI
STRY.fork0.chnk0: Indexing genome...
2018-04-28 20:04:49 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building transcriptome...
2018-04-28 20:04:58 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building kmer index...
2018-04-28 20:05:25 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Mapping reads...
2018-04-28 20:05:28 [runtime] (chunks_complete) ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY
2018-04-28 20:05:28 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.join
2018-04-28 20:05:31 [runtime] (join_complete)   ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY

2018-04-28 20:05:34 [runtime] (ready)           ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (split_complete)  ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS.fork0.chnk0.main
2018-04-28 20:05:40 [runtime] (failed)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS

>[error] The chemistry was unable to be automatically determined. This can happen if not enough reads originate from the given reference. Please verify your choice of reference or explicitly specify the chemistry via the --chemistry argument

What is going wrong here? How could I correct it?

You presumably aligned against the wrong reference transcriptome.

It is very weird: the thing that the first time I ran it, it worked. Then, because we got bad results after running cellranger count ('bad' means biologically not what we want to see. It would mean that the biology was done in the wrong way), we transferred the raw data again, generated fastq files again and somehow now I am having this issue for the first replica, but not for the second one which ran successfully. I am wondering where in this chain of actions I did a mistake. Can it be that the data got corrupted while I was transferring it, however, fastq generation worked still?

The issue was that I supplied wrong indexes during fastq files generation, however, they were generated successfully, and only the next step failed with not that informative error.


text:
Convert rs ID of one hg build to rs IDs of another build
I have a list of dbSNP rsIDs for GRCh37 and I want to convert them to the equivalent IDs in GRCh38. This is using the most recent dbSNP build (150 as of the time of this post). Is there any ID mapping available? If not, are there any tools I can use?

Did you try the liftOver tool of UCSC? You'll need a BED file with your SNPs coordinates for this.
You can also do this in R with rtracklayer.
library(rtracklayer)
?liftOver

There's no need to change the coordinates, the question is regarding the IDs themselves.

What I have is a  list of dbSNP IDs. I have downloaded source("http://bioconductor.org/workflows.R") workflowInstall("liftOver"), but not sure if ths will accept just a list of IDs.

No it won't accept only rsIDs, so you'll need to get the coordinates first with e.g., biomaRt. Sorry I understood that you had coordinates of the SNPs.

Thanks @b.nota I was able to get coordinates from biomaRT.

You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's chromosome|position|ID|weight, where position is sometimes empty, weight is typically 1, and ID is the rsID without the rs prefix.
Also, some rsIDs have a different strand. This is less likely to be an issue for you, but it's there in case you need it.

Thanks so much @Devon Ryan for all the useful information.

See also CrossMap:
http://crossmap.sourceforge.net/
20, 21, 22, 23, 24, 30

These are only the base positions that have changed between GRCh37 and GRCh38, not the IDs. For example if a study (e.g. found via pubmed) refers to a particular rsID, the genome build does not matter.
rsIDs do not depend on a reference genome. They point to a specified locus regardless of the differences in genomic assemblies. As a result, you cannot convert rsIDs of one genome build to rsIDs of another build. And there is no such mapping tool available, even if they pretend to be. rsIDs are assigned by dbSNP maintainers, and are stable. The numbers may only change on merge, where multiple numbers were assigned at the same genomic location, and the higher number gets merged into a lower number by the dbSNP maintainers, unless there is a very well-published, notable rsID mentioned in the press, in this case that merged number is considered the mentioned one, not the smaller one.

RefSNP number is the stable accession regardless of the differences in genomic assemblies.

For more information go to the “RefSNP Number Stability” section at https://www.ncbi.nlm.nih.gov/books/NBK21088/
The rsIDs may also be deleted in some cases see the following links on the explanation:

https://www.ncbi.nlm.nih.gov/books/NBK44467/#Build.i_found_rs123456_in_the_current_bu
https://www.ncbi.nlm.nih.gov/books/NBK44468/#Build.the__rscurrent_field_in_rsmergearc

Thus, specifying the location by rsID is stable, and can be only affected by an eventual merge or delete from the dbSNP database, if any, not by a new genome build. You can only “update” the rsIDs from time to time to make sure they all are the main rsID, not those larger numbers already merged to a smaller one. It should be noted that the base positions (offsets) of the same genomic location may change between the builds of the reference genome. And they did surely change between GRCh37 and GRCh38.
If your list is relatively small, you can get most recent rsIDs for your chromosome and offset via dbSNP API. It is available at https://www.ncbi.nlm.nih.gov/snp/docs/eutils_help/
The query for the API is the same as for the “advanced search” of dbSNP at https://www.ncbi.nlm.nih.gov/snp/advanced
If you specify data in a “Base Position” line, it will look for GRCh38 base position. If you need to specify the base position for a GRCh37, use “Base Position Previous” line. You may also need to add NOT MergedRS[All Fields] to exclude rsIDs that have already been merged.
For example, to get the rsID for a chromosome 1 GRCh37 base position 249222527, the API query URL is the following https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=1[Chromosome]+AND+249222527[Base+Position+Previous]+NOT+MergedRS[All+Fields]
You can do queries from a script if the number of requests is relatively small. You may find a perl script, for example, that does the job.
Even though the rsID numbers can be merged or deleted, and the changes are reflected in a new dbSNP database build, given the fact that the changes in numbers do not depend on the genome build, is better to always rely to the rsID number from the latest dbSNP database build.
The changes in rsID numbers are changed due to merges or deletions between dbSNP builds relatively infrequently, given the total number of the IDs in the database.


text:
bedtools feature out of bounds
I'm trying to count the GC in regions leading up to polyA clusters. I found the bedtools command nuc, one of the outputs of which is the GC content.
I took the mouse polyA cluster database from here, and the mm10 primary assembly soft masked from here. I imported the polyA clusters bed file to R. I then got a list of all the polyA clusters in the genes I'm interested in. For positive strands on polyA I took the start position, took 100 away from it and made that the start of the region of interest, then made the start of the polyA signal the end. For things on the negative strand I added 100 to the end of the polyA signal. I then exported the results to a bed file format. The code I used to do this is below.
When I try to run:
bedtools nuc -fi Mus_musculus.GRCm38.dna_sm.primary_assembly.fa -bed PolyA.bed
Every signal returns some variant of: Feature (chr10:21,080,315-21,080,415) beyond the length of chr10 size (0 bp).  Skipping.
I can however look at several of these regions in the UCSC genome browser without a problem, at least I don't think I had a problem.
Can anyone point me to what's going wrong or to another way I can solve this.
GetRanges <- function(FilteredClusters){
  for (x in 1:nrow(FilteredClusters)){
    if(FilteredClusters$Strand[x] == '+'){
  FilteredClusters$RegionEnd[x] <- FilteredClusters$ClusterStart[x]
  FilteredClusters$RegionStart[x] <- FilteredClusters$RegionEnd[x] - 100
}
else if(FilteredClusters$Strand[x] == '-')
    {
      FilteredClusters$RegionStart[x] <- FilteredClusters$ClusterEnd[x]
      FilteredClusters$RegionEnd[x] <- FilteredClusters$RegionStart[x] + 100
    }
  }
    ClusterCut <- select(FilteredClusters, 'Chr', 'RegionStart', 'RegionEnd')

    return(ClusterCut)
}

There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.

Well...that is depressingly simple. Thanks for the tip. It's worked just fine.


text:
Create galaxy local toolshed using Planemo
I tried to create galaxy local toolshed using planemo. I used command as follows,

create virtual env
virtualenv .venv; .venv/bin/activate

Install planemo
pip install planemo

To configuring a shed account
planemo config_init

Create a repository
planemo shed_init --name=seqtk_seq \
                  --owner=planemo \
                  --description=seqtk_seq \
                  --long_description="Tool that converts FASTQ
                                      to FASTA files using seqtk" \
                  --category="Fastq Manipulation"
planemo shed_create --shed_target local

This produces the following error:
Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.
(.venv) varuna@varuna-Inspiron-15-3567:~/git/galaxy_stable/planemo

Hi, could you paste the text of the error instead of screenshots? It might make it easier to read and help you.

Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.

Could not update local
HTTPConnectionPool(host='localhost', port=9009): Max retries exceeded with url: /api/repositories?key=dc5c7cd5cf56a2f44aacf755109e63b9 (Caused by NewConnectionError('<urllib3.connection.HTTPConnection object at 0x7f55542b7310>: Failed to establish a new connection: [Errno 111] Connection refused',)), 0 attempts left:

Try providing a project name planemo shed_create --shed_target local someProjectName.

I created galaxy local toolshed editing the xml files inside galaxy folder and  ran it using 9009 port. After run toolshed, I tried above mentioned command with project name.
planemo shed_create --shed_target local /home/varuna/git/galaxy_stable
Got following error.
"description required for automatic creation or update of shed metadata"

--shed_target local will only work if you're running a Galaxy toolshed on your local computer, since it's reading the .shed.yaml file you made with planemo shed_init and trying to create a repository with that name or a local (in this case) server. Planemo itself won't start a Galaxy toolshed, since it's fairly unusual for people to want to do testing involving creating local repositories. What you probably want to do instead is use planemo test instead. Once the planemo testing passes then use the various shed* commands to upload your wrapper to the toolshed or testtoolshed.


text:
How to read .bwt file .sa file?
In BWA project in bwt.c file I found bwt_dump_bwt method and bwt_dump_sa method. I want to use them to read .bwt file and .sa file. I wrote the following program:
int main(void)
{
  char *tmp, *prefix;
    bwt_t *bwt;
    prefix = bwa_idx_infer_prefix("hg19bwaidx");

    tmp = calloc(strlen(prefix) + 5, 1);
    strcat(strcpy(tmp, prefix), ".bwt");
    bwt = bwt_restore_bwt(tmp);

  return 0;
}

I got the following error:
Undefined symbols for architecture x86_64:
  "_bwa_idx_infer_prefix", referenced from:
      _main in readbwt-3022ec.o
  "_bwt_restore_bwt", referenced from:
      _main in readbwt-3022ec.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Any help ?
Thanks.

In the bwa source code, file bwt.{c,h} implement basic operations on .bwt and .sa files. They depend on utils.{c,h}, but you can create a fake utils.h to cut this dependency as is shown below. Then you only need three files: kvec.h, bwt.h and bwt.c to access .bwt and .sa.
The following example demonstrates how to find super-maximal exact matches (SMEMs). Note that internally bwa concatenates all database sequences into one long string. If you use bwt.{c,h} only, you get coordinates on this long string but not on each chromosome. You have to call APIs in bntseq.{c,h}, which are fairly independent, too. You can go on to implement a mini bwa-mem with this example, but you will need all the bwa-mem source code, except the main function.
Compiling and running
# compile the example
git clone https://github.com/lh3/bwa
mkdir -p tiny-smem && cd tiny-smem
cp ../bwa/{bwt.c,bwt.h,kvec.h} .              # you need at least these three files
### then create utils.h and test.c shown below ###
gcc -g -O2 -Wall bwt.c test.c -o tiny-smem    # compile tiny-smem

# create a bwa index if you don't have it. We use chr20 from GRCh38 as an example:
(cd ../bwa && make)                           # compile bwa for indexing
../bwa/bwa index -p chr20 path/to/human/chr20.fa

# run the example
./tiny-smem chr20.bwt chr20.sa TTTTCTTCTTTTGTAATATAAGCATTgGTTGTTCTAAATTTCTCTCTTGGCACTGCTTTA

The last command line above outputs:
0       0       26      1       -125734
0       14      27      4       -3657204        +16021887       +11802682       +54692983
0       15      29      1       +8344885
0       17      31      1       -39545379
0       19      37      1       +55060482
0       25      39      2       +1198672        -26409692
0       27      60      1       -125700

where column 2 and 3 give the coordinate of SMEM on the query, col4 the number of exact hits and the following columns are the list of hits in the concatenated coordinate.
File: utils.h
#ifndef FAKE_UTILS_H
#define FAKE_UTILS_H

#define xassert(cond, msg) assert(cond)
#define xopen(fn, mode) fopen((fn), (mode))
#define err_fclose(fp) fclose(fp)
#define err_fread_noeof(ptr, size, n, fp) fread((ptr), (size), (n), (fp))
#define err_fwrite(ptr, size, n, fp) fwrite((ptr), (size), (n), (fp))
#define err_fflush(fp) fflush(fp)
#define err_fseek(fp, off, whence) fseek((fp), (off), (whence))
#define err_ftell(fp) ftell(fp)

#endif // ~FAKE_UTILS_H

File: test.c
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include "bwt.h"

int main(int argc, char *argv[]) {
  bwt_t *bwt;
  int j, max_intv = 100; // ignore SA intervals larger than this value
  bwtintv_v r = {0,0,0};

  if (argc < 4) { fprintf(stderr, "Usage: ./test <.bwt> <.sa> <seq>\n"); return 1; }
  if ((bwt = bwt_restore_bwt(argv[1])) == NULL) abort(); // read .bwt
  bwt_restore_sa(argv[2], bwt);                          // read .sa
  for (j = 3; j < argc; ++j) {
    int i, k, x, len = strlen(argv[j]);
    uint8_t *seq = (uint8_t*)malloc(len);
    for (i = 0; i < len; ++i) { // convert A/C/G/T/N => 0/1/2/3/4
      int c = tolower(argv[j][i]);
      if (c == 'a') seq[i] = 0;
      else if (c == 'c') seq[i] = 1;
      else if (c == 'g') seq[i] = 2;
      else if (c == 't') seq[i] = 3;
      else seq[i] = 4;
    }
    x = 0, r.n = 0;
    while (x < len) {
      x = bwt_smem1(bwt, len, seq, x, 1, &r, 0); // find SMEMs
      for (k = 0; k < r.n; ++k) { // traverse each SMEM
        bwtintv_t *p = &r.a[k]; // p->info encodes start and end on the query
        printf("%d\t%d\t%d\t%d", j-3, (int)(p->info>>32), (int)p->info, (int)p->x[2]);
        if (r.a[k].x[2] <= max_intv) { // don't retrieve coordinates for huge intervals
          int l_mem = (int)p->info - (int)(p->info>>32);
          for (i = 0; i < p->x[2]; ++i) { // traverse each hit
            uint64_t y = bwt_sa(bwt, p->x[0] + i); // get suffix array coordinate
            if (y < bwt->seq_len>>1) printf("\t+%ld", (long)y);
            else printf("\t-%ld", (long)(bwt->seq_len - y - l_mem));
          }
        }
        putchar('\n');
      }
    }
    free(seq);
  }
  free(r.a);
  bwt_destroy(bwt);
  return 0;
}

You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c.
Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.

Thanks. But there is no main() in bwa.c by the way I want to do something ppl did before. like this call to bwt_dump_bwt in this link: https://github.com/nsoranzo/bwa/blob/master/bwtindex.c The problem is I'm new to c programming. That why I need a help and still I can't understand how to do it.

There's a main in a different file then. You linked to a fork of the bwa repository, which apparently used to include that file.

I ran clang -o main bwa.o main.c -lz -lpthread -lm -lrt and I got: ld: library not found for -lrt and when I removed  -lrt I got many other errors. Any idea? I'm a Mac user maybe this is the problem.

You'll have to figure out the equivalent library (or install it) on OS X, or whatever system you're using.

I ran BWA successfully on my system. But I can't read .bwt with the program that I wrote.


text:
Best data repository to publish a large 'plotly' table containing all annotations on a transcriptome?
I am not sure whether this is the best SE community for this question.
I want to publish a transcriptome paper, along with interactive materials enabling readers to peruse the data behind the discussion. The R package 'plotly' enables rendering online-publishable tables that can be filtered & searched for specific entries.
I would like to publish such a table as supplementary to a paper (e.g. preprint version) however usually data repositories allow for uploading files and text.
Would anyone here know an official data repository with a format such that would enable publishing online an interactive table?
A good example of such a table is given here.
Thanks in advance.

This might be better suited to [Academia] where people from more fields could know about a place for such things.

@Llopis I will also ask there, thanks! My guess is that they'll think it is mostly technical and redirect me elsewhere. I think it is in fact a quite specific question on a situation not many have faced yet.

You may try Code Ocean. Currently, it is not that easy to use IMHO, but one of their tech leads promises to substantially improve it towards the year end.

I have just looked into Code Ocean as you suggested, and though it looks promising it is indeed not very user-friendly at the moment... Plus, looks like one needs a subscription to run interactive content.

I have seen some (1 or more) articles with interactive data displays in the F1000Research journal I don't know recall now if they had interactive tables or plots but they used plotly for sure.

Yes, I heard that F1000 encourages plotly content. I have just sent them an email. However that would be a journal option for a final publication (if I can pay their fees). Ideally, right now I'd greatly benefit from an external repository that could also be used with a preprint server like Bioarxiv. Unfortunately bioarxiv seems to only take .pdf output.

Yes, I realized later that your data seems "huge" so a publisher doesn't seem a good choice to store the data

Mind that by "huge" we mean a table is very large in numbers of rows. In reality the actual data behind the table are .txt files adding up to less than 20 Mb. I could in principle publish such table as a .pdf file, however nobody would like it and read it.

The example looks a lot like a shiny application. Your best bet then is one of the following:

Ask your IT department to setup a shiny server for you. They can use the free version.
Post your shiny app on shinyapps.io. Depending on how popular the page becomes, you might have to pay for an account.

The free version of shiny works for most cases (it's what we use), the only problem with it is that it doesn't scale (but if you care about that then you need to either pay for hosting on a pro server or buy a shiny pro license).

If you have problems scaling shiny apps you might find this post interesting.

Only now I see the issue you mention about popularity. On a free account I only get 24 mins of server availability per month. Given my scope of use, this is way too short !..


text:
bioawk removed part of FASTQ header
I used
bioawk -cfastx 'length($seq) > 1 {print "@"$name"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

in order to keep a particular read length, but this command shortened the header  from
@A00199:161:HF3JLDMXX:1:1101:5882:1063 2:N:0:CTGTAGCC+GAACAGGC

to
@A00199:161:HF3JLDMXX:1:1101:5882:1063

How is it possible to fix the bioawk command to remain the full command?

bioawk -cfastx 'length($seq) > 1 {print "@"$name $comment"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

Note the addition of $comment, which is what follows the name on the first line of each entry.

Space was missing and I had to change it to {print "@"$name " " $comment"\n"$seq"\n+\n"$qual}


text:
STAR Indexing Diference for Small and Large Genome File Output
I have a quick doubt on the output of the Genome Indexing, I have used the STAR program along with genome .fasta file and GFF file.
Genome size is 3GB, here is the file output
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
genomeParameters.txt

I have another small Genome 60MB in size, I did the genome indexing, here is the file output
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
exonGeTrInfo.tab
exonInfo.tab
geneInfo.tab
Genome
genomeParameters.txt
SA
SAindex
sjdbInfo.txt
sjdbList.fromGTF.out.tab
sjdbList.out.tab
transcriptInfo.tab

My point here is that, why I got the extra information for my small genome size, but I didn't get the same for the big size genome. I do apply the same procedure for the both.
here is the below information. Only difference I made for the large Genome size is (--sjdbOverhang 99 --genomeChrBinNbits 15) to reduce the memory, but the rest of things are same for small genome.
#!/bin/bash
NUM_THREADS=12
mkdir DB
STAR --runMode genomeGenerate --genomeDir DB \
    --runThreadN $NUM_THREADS \
    --genomeFastaFiles XL9_2.fa \
    --sjdbGTFfile XENLA_Frog.gtf+gff3 \
    --sjdbOverhang 99 \
    --genomeChrBinNbits 15

Could anyone give an idea, why there is different, I am new to this field, so I am wondering about the difference in this.
Thanks in advance.
Cheer
San

The files missing in the case of the larger genome are the actual index itself - it doesn't look like STAR completed. Do you have STDOUT/STDERR from the runs?

Here is my below syntax, I use it for It                                                     #!/bin/bash
module load star/2.5.0c
NUM_THREADS=22
mkdir Sangv92db
STAR --runMode genomeGenerate --genomeDir Sangv92db
--runThreadN $NUM_THREADS
--genomeFastaFiles /ifs/home/skm9/san/xlaevis/XL9_2.fa
--sjdbGTFfile /ifs/home/skm9/san/xlaevis/XENLA_9.2_Xenbase.gff3
--sjdbOverhang 99
--genomeChrBinNbits 15,                                                                                It run and gene the files out same like the above I mentioned.

here is the memory allocation I used for it, still not able to complete the indexing, ... qsub -pe threaded 4-20 -hard -l mem_free=400G -l h_vmem=560G -l mem_token=80G RunSTARIndex2.sh...

You didn't get extra information with the small genome, the index for the larger one just wasn't made. My guess is that you ran out of RAM. For mammalian sized genomes ensure you have at least 30GB of RAM, likely more like 60GB. Also, make sure to look for error messages.

qsub -pe threaded 4 -hard -l mem_free=60G -l h_vmem=80G -l mem_token=60G RunSTARIndex.sh...do you think this much is enough, although one node has 128GB....hope this time it will....me running now...Thanks a million in advance...

Thanks for the information, I do run the same but it seems, its terminating and giving the same output files. ..this is what I tried in the last ...qsub -pe threaded 4-20 -hard -l mem_free=100G -l h_vmem=120G -l mem_token=60G RunSTARIndex2.sh....Each Node -128 RAM, is their other possible way to Index it?

Some genomes require gobs of memory. If you point me to the files you're trying to index I can try and do that for you.

Thanks a million Devon Ryan, I tried all possible way to increase memory like ....qsub -pe threaded 4-20 -hard -l mem_free=400G -l h_vmem=560G -l mem_token=80G RunSTARIndex2.sh, but still fail, I don't see any error message, its seems like it terminate on the half way...Could you let me know your email, My be I will email the same to you.

Genome Indexing for Xenopus leaves, v9.2, here is the link for the database http://www.xenbase.org/other/static/ftpDatafiles.jsp, Please let me know, if you can  do the indexing for it, Let me know, how I can transfer the index file too, thanks a million in advance ...


text:
how to export the monocle data into Seurat
I have tried to export the monocle data into Seurat and got the error
"Error in Seurat::CreateSeuratObject(raw.data = data, normalization.method = "LogNormalize",  :
  unused arguments (raw.data = data, normalization.method = "LogNormalize", do.scale = TRUE, do.center = TRUE, is.expr = monocle_cds@lowerDetectionLimit)
Calls: exportCDS
In addition: Warning message:
In if (export_to == "Seurat") { :
  the condition has length > 1 and only the first element will be used"

Here is my code seurat_SAN_cds_subset <- exportCDS(SAN_cds_subset)
Is anyone familiar with this?

There's no default for the export_to option, so you just need to specify export_to="Seurat". Note that monocle2 doesn't seem to be compatible with the newest Seurat (3.1) releases.

It is possible to update Seurat v2 objects to Seurat v3 objects via UpdateSeuratObject().

I don’t know off hand, maybe give it a whirl and see.


text:
SC-RNA seq percent.mito
Is that possible that the percent.mito level is zero? I have done vlnplot and it shows that its percent.mito level is zero. Does that mean my  mitochondrial gene didn't detect?

It is violin plot. And the command is:

Sorry hua, but I find it difficult to understand. Could you clarify what is "percen.mito" ? It seems like it should be a function or a property of your data but I'm not sure. Also, could you explain what is "VlnPlot" ? It seems like a violin plot, but I'm not sure, perhaps the image didn't get attached (upload just a link and we will add it).

Your genome likely doesn't title mitochondrial genes MT-something, so you'll need to change how mito.genes is created.

Is "percent.mito level" the percent of reads mapping to genes annotated as mitochondrial? Are the mitochondrial genes annotated? Note you have [edit] button you can use to improve your question.

You should put "Seurat" in your question, and add a Seurat tag.

I've just come across a similar issue. It seems the mouse and human mitochondrial genes use mixed capitalizations. See the following PR for details:
https://github.com/Hoohm/dropSeqPipe/pull/43
Bottom line is, if you ignore the case in the grep pattern it might work (assuming you use human/mouse):
mito.genes <- grep("^mt-", rownames(shox2crepos_mtmg@data), value=TRUE, ignore.case=TRUE)

Thank you so much, I have solve this problem after your instruction.

The simplest answer is that the bit of code you copied from the Seurat vignette is not working on the gene names you have.  Look at the code.  If none of your gene names begin with "MT-", none will be identified as mitochondrial.


text:
Assembly by stringtie
I run this cmd
./stringtie G1_sorted.bam -B -o G1.gtf -G Triticum_aestivum.IWGSC.42.gtf -p 4  -C G1.refs.gtf -A G1.abund.tab

Error is:
-WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
Please make sure the -G annotation file uses the same naming convention for the genome sequences.

How can I resolve it?
Update:
I have downloaded both (GTF and Refernce Genome)files from Ensemble. Can you confirm me that whether theu have same chromosomes naming? So that i can perform my mapping  task of GTF File with Reference Genome which usually i remain fail to do due to difference in chromosomes naming.                                                           Here The Reference genome file is shown as:
 >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC

Where as The annotation file (gtf file ) is as
    #!genome-build IWGSC
#!genome-version IWGSC
#!genome-date 2018-07
#!genome-build-accession GCA_900519105.1
3B      IWGSC   gene    212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; gene_source "IWGSC"; gene_biotype "protein_coding";
3B      IWGSC   transcript      212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; transcript_id "TraesCS3B02G000100.1"; gene_source "IWGSC"; gene_biotype "protein_coding"; transcript_source "IWGSC"; transcript_biotype "protein_coding";

Can it be confirmend that these both files are same ? If not same with chromosome naming convention , then please suggest me how can I change it ?

This looks like your chromosome names in the BAM don't match the GTFs you're using. Check that first.

the reference genome download from the NCBI.and gtf  file from ensembl. your mean is that i have to use the gtf file from the NCBI, then this error will be removed?

Welcome to the site @fatima. Yes ENSEMBL uses genes identifiers starting with ENSGXXXXX.XX while NCBI uses gene names with just numbers and they don't map to the same positions. They aren't exacts equivalent

hi IIrs, i have use the reference genome and gtf file from ensembl.encounter same problem.how i can mofify the file with same name?

hi @fatima (to send me a notification use @ followed by my username) Could you please explain how do you want to modify these files to use the same name? rereading again your question I find it a bit unclear what do you want to achieve and what do you expect us to do for you.

This seems to be a common problem. As @llrs pointed out, Ensembl, UCSC(genome browser and associated data) and NCBI records are not directly interchangeable. The ambiguity refers to annotation as well as the reference genome. The reference genome is identical, except for the chromosome naming convention, meaning a reference from Ensembl might not contain CHR prefix, whereas it will be present in one downloaded from UCSC.
To prevent such conflicts in downstream tools, it is good to stick to the reference and annotation from the same repository throughout the analysis and map the genes to gene names or whatever towards the end of the investigation.
You can download all the associated data for T. aestivum which I beleive is the organism you want to study from Ensembl plants FTP download.

Hi Siddharth, Thank you, Now I have downloaded both files from Ensembl. Can you please tell me these are with the same naming?                                                                                I have open the reference file.The Reference genome file is as                                                >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC
The annotation file (gtf file ) is as


text:
GEOquery errors due to list input
I want to extract the expression matrix from the GSE object in GEOquery.
library(GEOquery)

gse1 <- getGEO("GSE97168",GSEMatrix=TRUE)
datExpr = exprs (gse1[[1]])
datExpr

GSM2553017 GSM2553018 GSM2553019 GSM2553020 GSM2553021 GSM2553022 GSM2553023 GSM2553024 GSM2553025
     GSM2553026

Table(GSMList(gse1)[[1]])[1:5,]

Table(GSMList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GSMList’ for signature ‘"list"’

Table(GPLList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GPLList’ for signature ‘"list"’

Meta(gse1)$description
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Meta’ for signature ‘"list"’

Columns(gse1)
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Columns’ for signature ‘"list"’

Why it's not working with R version 3.5.1?

Welcome to the site. Are you sure this should be working? In the vignette Table is only used to "Look at data associated with the GSM" not to retrieve the data expression data or information about a specific GSM

Use GSEMatrix=FALSE, the information you want appears to be in the soft file, rather than the matrix. Yes, this is moderately annoying since it means you need two objects to describe this particular dataset.


text:
Geneland R + Fortran error?
I am helping an undergrad with their project and their thesis advisor has requested they use GENELAND for their analysis. I haven't worked on this sort of analysis, so all I have done is install it in their R, hand over the manual, and help troubleshoot any errors they can't resolve.
We've gotten the analysis through the MCMC stage and are getting an error with the PostProcessingChain command.
Their MCMC command is:
MCMC(geno.hap=geno1, spatial=FALSE, nb.nuclei.max=3,
    varnpop=TRUE, freq.model="Uncorrelated",
    filter.null.alleles=FALSE, nit=100000,
    thinning=100, npopmin=1, npopinit=3,
    npopmax=10, rate.max=3, delta.coord=0,
    path.mcmc="/Users/myname/localdoc/GENELAND")

Their PostProcessChain command is:
PostProcessChain(coordinates=gene1coord.utm,
    path.mcmc ="/Users/myname/localdoc/GENELAND/",
    nxdom=100, nydom=100, burnin=200)

This is exactly the same command as in the manual, and gene1coord.utm is an R object with the coordinates transformed into UTM coordinates. The error we receive is:

([1] "Calling Fortran function postprocesschain2"

At line 12334 of file main.f (unit = 10)

Fortran runtime error: File already opened in another unit)

I have tried fix(PostProcessChain) to see if there's something I can adjust and this is what I see for the section that shows up:

So far, here is the list of things I have tried:

Reinstall X11
Reinstall gfortran
make sure that R is 3.4.x as the GENELAND package suggests
Tried different parameters on MCMC just in case it was user error
Tried different parameters on PostProcessChain, but they're the same as the manual really
Tried different kinds of coordinates
Tried a spatial analysis instead of a non-spatial
Googled the error
Restarted the computer

Alternatively, is there another program that does similar analysis to geneland? I'm not sure if the advisor recommended it because it's the only thing that does what it does or because that's just what they used back in the day.
Edit: There were 3 blank files in the folder titled %FO, P, and x. When I deleted them the "file open" error went away but it still hangs up on [ 1 ] Calling Fortran function postprocesschain2 and has a fatal error that crashes R.

I don't suppose the file is currently being accessed by anything else, right? Just to be 100% sure, have you tried rebooting (in case something weird happened and you have a stale open file handle somewhere)? Also, please avoid posting images of text. That makes it considerably harder for us to copy the error and look for it.

@terdon Thanks, no one/nothing else could be accessing the files but R. There isn't even a file called "main.f" or any files with 12000 lines in the folder the MCMC output goes to, and PostProcessChain accesses.
I tried a restart and it changed nothing.


text:
How I normalize these two sets of data
I have average log fold change for a cluster of cells versus another cluster of cells by Seurat like below
         avg_logFC
TAGLN   -2.507833541
TXNRD1  1.093672244
MYL9    -1.579165463
ACTB    -1.482386539
ACTA2   -2.379757433
SQSTM1  1.227562767
ACTG2   -2.144757086

I have also mean log2 ratio of a set of proteins in treatment versus control like below (I converted protein IDs to gene symbol)
mean log2ratio  Gene Symbol
-0.3    CYFIP1
0.4     GNPTG
-0.3    WWC2
-0.5    TRPV4
-0.4    CHRM2
0.3     MFSD11
-0.3    ANXA11
0.3     CENPE
-0.4    RPL22L1
0.3     SAMD1

How I can normalize these two data or make them compatible for plotting a heatmap, pathway analysis, etc?

You cannot normalise this further because it is already transformed and standardied "mean log2 ratio".  Sauret uses PCA and tSNE, where normalisation is not needed. You have previously demontrated advanced skills in machine learning, for example using neural networks, so I am very surprised by this question

There's no need to normalize them, you're not comparing them. Just use them as they are.

Thank you, even if one comes from single cell RNA-seq and another comes from proteomics? The ultimate goal would be visualizing their relationship

The sources doesn’t really matter. You can scale them a bit if it helps for visualization, but that’s as much as would be useful.

Sorry to coming back again. If I want to combine these two data sets; for instance, I would have gene A common between two data sets. So, If I want to visualize both data in combination not comparison, how I deal with two values from different experiments for this overlap genes? Can I simply take average of values?

That's up to you, there's no single answer to such questions.


text:
Changing the colour of each cell in tSNE plot
I have plotted a tSNE plot of my 1643 cells from 9 time points by seurat like below as 9 clusters. But, you know I should not expected each cluster of cells contains only cells from one distinct time point. For instance, cluster 2 includes cells from time point 16, 14 and even few cells from time point 12. My question is, how I can change the colour of cells here so cells from each time point has a distinct colour does not matter some cells are being clustered somewhere else because this is unavoidable. I have named my cells already so I know cells from each time point but I want to colour cells to nine colours based on the time point they have been sequenced not the cluster they are being grouped by seurat.
I first created a file in one column is name of my cells and in second colour the corresponding time point
Like so
    time
s0.1    h0
s0.2    h0
s0.3    h0
s0.4    h0
s0.5    h0
s0.6    h0
s0.7    h0
s0.8    h0
s0.9    h0
s0.10   h0
s0.11   h0
s0.12   h0
s0.13   h0
s0.14   h0
s0.15   h0
s0.16   h0
s0.17   h0
s0.18   h0
s0.19   h0
s0.20   h0
s0.21   h0
s0.22   h0
s0.23   h0
s0.24   h0
s0.25   h0

I read my file
time <- read.table("time.txt",header = TRUE)
seurat@meta.data$time <- time
TSNEPlot(object = seurat, group.by = "time")

Sorry in my code colour of tSNE changed also the position of cells changed that I don't want that. I want to change the color of cells but they are still in the same cluster
says that
> TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
  'to' must be a finite number

when I changed names(identity)=sampleInfo$cellName to names(identity)=src$cellName
that worked but again the position of cells changed that I don't want that :(

If you can extract the underlying data (as this is essentially an xy scatter plot) you can then share it here.  Your question then essentially becomes a base R / ggplot2 problem which should enable others to help you better than the screen shot alone.

From an experimental interpretation point of view if your cells are failing to cluster by time point then clearly there is something else going on with your data, again more information would be useful, additionally if you want help with plotting we'd need to know which time point each data point is from.  Thus some exemplar data is required for us to be able to help.

Please  no negative point :( If I lose my reputation I can not ask question anymore

I found the solution   if this is my time table data cc <- AddMetaData(object = n, metadata = time)
cc <- SetAllIdent(object = cc, id = "time")
TSNEPlot(object = cc)
This only colours the cells based on the time point they are

Seurat puts the label in the tSNE plot according to the @ident slot of the Seurat object.
By default, if you do the tSNE without computing the clusters and you have the correct metadata in the object, the labels should be pointing to your timepoints not to the clusters.
If it is not the case, what you have to do its to change this slot to reflect to which timepoint belong each cell you are plotting.
For instance if your Seurat object is called src:
identity=as.factor(sampleInfo$timepoints)
names(identity)=sampleInfo$cellName
src@ident=identity

And then plot the tSNE. Something like this should work after adapting it to your metadata.

Sorry what is sampleInfo?

What you call time in your example. cellName could be your s0.1, etc and timepoints your h0, etc

Sorry says that  > TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
'to' must be a finite number

Sorry I can nit understand. I have the name of my cells and time points as I showed in my question, so sampleInfo is name of my time contains cell names and time points?

Yes, it is supposed to contain similar information. I put this name randomly. In your case it would be your time data frame

Sorry but says that    > TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
'to' must be a finite number