text:
What is "Possible polymerase run-on fragment"?
I used gffcompare recently to compare a newly assembled transcriptome and a published transcriptome. The class code "p" transcripts are "possible polymerase run on".

What does that mean? I tried very hard but still cannot find the answer.

The general idea behind a run-on fragment is that it's background noise. This derives from an open area of the genome that's next to an area that's actually being transcribed. Thus, all of the machinery for transcription is already present, so if the region adopts a bit more of an open conformation then that machinery can sporadically bind and start transcribing. You can get those either up- or down-stream of genes. A lot of things like this are just background transcriptional noise.

Thanks, Devon!!! Is there a paper, for example, a review or research paper related to this topic that I can read more (and cite in the future)?

Not that I know of, though perhaps @IanSudbery knows of one.

On most coding genes (with the exception of replication dependent histone genes), Transcripts (as opposed to transcription) are terminated by cleavage polyadenylation sites, where the growing transcript is cleaved by a protein complex and the poly-A tail. However, PolII continues to transcribe until it is terminated by torpedo termination (see here). Normally this extra RNA after the polyA site is degraded by Xrn2, but if you capture any particular snapshot of the RNA in the cell at any one time there will be a small amount that corresponds to this run on transcription.

Thanks for your help, Ian! The situation you mentioned seems could only occur downstream of the gene. Is that correct?

Yes, this should only happen downstream of a gene.


text:
What does an output of -I for all amino acids mean in a psi-blast pssm matrix file?
I have run psi-blast using the NR database, remotely, with one iteration, for several sequences to calculate an evolutionary profile (PSSM) for each of those sequences.
However, many of the PSSM files contain lines of -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I rather than numerical values.  What does it mean?
>     Last position-specific scoring matrix computed
>                 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
>         1 G  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>         2 N     4  -1  -2  -2   0  -1  -1   0  -2  -1  -1  -1  -1  -2  -1   1   0  -3  -2   0
>         3 K    -2  -1   4   4  -3   0   1  -1   0  -3  -4   0  -3  -3  -2   0  -1  -4  -3  -3
>         4 E     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>         5 K    -2  -2   1   6  -3   0   2  -1  -1  -3  -4  -1  -3  -3  -1   0  -1  -4  -3  -3
>         6 A    -1   0   0   2  -4   2   5  -2   0  -3  -3   1  -2  -3  -1   0  -1  -3  -2  -2
>         7 D    -2  -3  -3  -3  -2  -3  -3  -3  -1   0   0  -3   0   6  -4  -2  -2   1   3  -1
>         8 R     0  -2   0  -1  -3  -2  -2   6  -2  -4  -4  -2  -3  -3  -2   0  -2  -2  -3  -3
>         9 Q    -2   0   1  -1  -3   0   0  -2   8  -3  -3  -1  -2  -1  -2  -1  -2  -2   2  -3
>        10 K    -1  -3  -3  -3  -1  -3  -3  -4  -3   4   2  -3   1   0  -3  -2  -1  -3  -1   3
>        11 V    -1   2   0  -1  -3   1   1  -2  -1  -3  -2   5  -1  -3  -1   0  -1  -3  -2  -2
>        12 V    -1  -2  -3  -4  -1  -2  -3  -4  -3   2   4  -2   2   0  -3  -2  -1  -2  -1   1
>        13 S    -1  -1  -2  -3  -1   0  -2  -3  -2   1   2  -1   5   0  -2  -1  -1  -1  -1   1
>        14 D    -2   0   6   1  -3   0   0   0   1  -3  -3   0  -2  -3  -2   1   0  -4  -2  -3
>        15 L    -1  -2  -2  -1  -3  -1  -1  -2  -2  -3  -3  -1  -2  -4   7  -1  -1  -4  -3  -2
>        16 V    -1   1   0   0  -3   5   2  -2   0  -3  -2   1   0  -3  -1   0  -1  -2  -1  -2
>        17 A    -1   5   0  -2  -3   1   0  -2   0  -3  -2   2  -1  -3  -2  -1  -1  -3  -2  -3
>        18 L     1  -1   1   0  -1   0   0   0  -1  -2  -2   0  -1  -2  -1   4   1  -3  -2  -2
>        19 E     0  -1   0  -1  -1  -1  -1  -2  -2  -1  -1  -1  -1  -2  -1   1   5  -2  -2   0
>        20 G     0  -3  -3  -3  -1  -2  -2  -3  -3   3   1  -2   1  -1  -2  -2   0  -3  -1   4
>        21 A    -3  -3  -4  -4  -2  -2  -3  -2  -2  -3  -2  -3  -1   1  -4  -3  -2  11   2  -3
>        22 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        23 D    -2  -2  -2  -3  -2  -1  -2  -3   2  -1  -1  -2  -1   3  -3  -2  -2   2   7  -1
>        24 M    -1   0   0   1  -3   4   4  -2   0  -3  -3   1  -1  -3  -1   0  -1  -2  -2  -2
>        25 Y     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
>        26 K    -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4  -4
>        27 L    -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
>        28 D    -1  -2  -3  -3  -1  -2  -3  -4  -3   3   3  -3   2   0  -3  -2  -1  -2  -1   2
>        29 N  -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I -I
>        30 S    -2   1  -1  -1  -2  -1  -2  -2  -1  -1  -1  -1  -2   0  -1   0   0  -1  -1  -3
>        31 R     4   0  -2   1   4  -1  -3   0  -3   0  -3  -4   4  -1   0  -1  -3  -1   0  -4
>        32 Y    -3  -1  -3  -4  -3  -3  -2  -3  -1  -3  -1  -1  -3  -3  -3  -1  -1  -1  -3  -2

I don't know but guess that either something has gone wrong or those somehow represent insertions/deletions (gaps)

can you try running with more iterations? Running with a single one makes little sense I think

@Chris_Rands Perhaps more iterations would help, however, remote blast only allows one iteration, unless requerying.  I could try running locally to see.  Strangely, for SWISSPROT database (tiny compared to NR) it doesn't appear to happen (the instances of -I).

It means -inf. It has to do with how the PSSM score is calculated.

That sounds very plausible, do you have a source/reference? Aren't pseudocounts normally added by default? will running >1 iteration change this?

Yes it appears to be minus infinity according to that wikipedia article.  Strange how it was so hard to find this information.  Thanks!  Now I just need to figure out what to do about it since I need to encode this data.

@Chris_Rands, if you have blast installed type in psiblast -help. You'll see that pseudocount is set by default to 0.

On NCBI PSI-BLAST web interface it says "Pseduocount parameter. If zero is specified, then the parameter is automatically determined through a minimum length description principle (PMID 19088134). A value of 30 is suggested in order to obtain the approximate behavior before the minimum length principle was implemented."

@Chris_Rands, did you read psiblast -help or not? That is something else than a web interface.

source: Stack Exchange

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

sort: votes


text:
Tumor sample and heterozygous SNVs: what is the genotype of the normal cells?
In the following picture, you can see that in normal cells (yellow) we have two copies of each gene, while in tumor cells (violet), we just have one copy due to a clonal mono-allelic deletion. I do not really care about the genotype of the two genes, I just need to know that they are heterozygous.

On the other hand, suppose we are talking about somatic heterozygous SNVs and therefore we do not consider deletions anymore. In this case I would like to know the genotype (AA, AB or BB) of tumor and normal cells. Based on my understanding, I expect the AF of the SNVs in the tumor cells to be roughly 0.5 (since I ask for heterozygous SNVs). But what about normal cells? What is their genotype? If I ask for heterozygous SNVs, does that imply that normal cells are homozygous for that alleles?

Normal cells also have SNV and mutations, in fact recognizing which mutations and SNV are the ones that cause the tumor is an active research area. But it is impossible to say which is their genotype, as you would need to sequence them to know it.

Your initial information/image and actual question seem pretty unrelated to me, but I'll do what I can.
If you have a heterozygous SNV in a malignant population of cells, it really depends on when the mutation occurred for you to know the genotype of wildtype cells. It could have been there since birth (a germline mutation passed down from the parents) or have been acquired throughout the person's life and clonally expanded (somatic mutation).
Unfortunately, you have no way to differentiate between these two types of SNVs just by sequencing the tumor cells. If you want to know true somatic mutations that are found only in the tumor cells, you need to sequence a matched normal control from the same patient (preferably from the same cell type). This is why most large cancer sequencing studies include not only tumor samples, but normal cells from the same patient - see TCGA or ICGC for examples.

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


text:
Normalize by counting CNV within a pool of VCF files
I’m (trying) using the GATK4 germline CNV calling pipeline. I successfully got 57 VCFs from my sample batch, called with segments (obtained by merging the contiguous intervals), like in a classic VCF :
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  2046745451-1006_S4
M       3288    CNV_M_3288_15907        N       <DEL>,<DUP>     .       .       END=15907       GT:CN:NP:QA:QS:QSE:QSS  2:5:9:17:6:21:21
1       69071   CNV_1_69071_70028       N       <DEL>,<DUP>     .       .       END=70028       GT:CN:NP:QA:QS:QSE:QSS  1:0:1:204:204:204:204

But I got way too much of those intervals, more than 10k. I would like to know if there is an existing tool which count the different segments (variants / intervals common by +/= 75% of their length) in one VCF and gives me the count of the different segments overlapped by segments in other sample in my batch. By counting the most redundant segment, I could determine which are background noise and maybe decrease the number of variants in my VCF by filtering.
Thank you.


text:
How to interpret a PSSM matrix from PSI-Blast - what do the numbers represent?
PSI-Blast on a sequence returns a PSSM matrix file containing two matrices, one to the left, one the the right, for all amino acids:
A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V

What is the difference between these two matrices, how can I interpret them?
Last position-specific scoring matrix computed, weighted observed percentages rounded down, information per position, and relative weight of gapless real matches to pseudocounts
            A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V   A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
    1 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    2 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    3 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    4 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    5 L    -2  -3  -4  -4  -1  -2  -3  -4  -3   2   4  -3   2   0  -3  -3  -1  -2  -1   1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    6 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    7 Y    -2  -2  -2  -4  -3  -2  -2  -4   2  -2  -1  -2  -1   3  -3  -2  -2   3   8  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    8 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
    9 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   10 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   11 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   12 F    -3  -3  -3  -4  -3  -4  -4  -4  -1   0   0  -4   0   7  -4  -3  -2   1   3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   13 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   14 V     0  -3  -3  -4  -1  -3  -3  -4  -4   3   1  -3   1  -1  -3  -2   0  -3  -1   4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   15 C     0  -4  -3  -4  10  -3  -4  -3  -3  -1  -1  -4  -2  -3  -3  -1  -1  -3  -3  -1    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   16 V     0  -3  -3  -4  -1  -3  -3  -4  -4   3   1  -3   1  -1  -3  -2   0  -3  -1   4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   17 G     0  -3   0  -2  -3  -2  -2   7  -2  -4  -4  -2  -3  -4  -2   0  -2  -3  -4  -4    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00
   18 R    -2   6  -1  -2  -4   1   0  -3   0  -3  -3   2  -2  -3  -2  -1  -1  -3  -2  -3    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0.00     0.00

Additionally, at the end of the file, there is additional information, what does it show exactly?
                      K         Lambda
Standard Ungapped    0.1628     0.3454
Standard Gapped      0.0410     0.2670
PSI Ungapped         0.1628     0.3454
PSI Gapped           0.0410     0.2670

How do you actually obtain the PSSM matrix?  Can you only get it from command line?  The user interface on NCBI only allows download of a .asn file.

@BrianWiley You can download the PSSM from the web search tools: make sure you select PSI Blast search and perform more than 1 iteration. The PSSM is only available after more than 1 iteration. Moreover, you can instead download and install NCBI BLAST+ suite from https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download ...once installed, with a protein database also downloaded, also available from that page (usually people use NR, swissport or uniref), then you can search with the psiblast program.  Note: If you want to use uniref you will usually have to generate it.

Thanks.  How much space is required for the NR database?  I was able to get this format of the profile using psiblast with command line like you indicated.  Just made my own small database.

If you read wiki carefully, you can see that PSSM is calculated in 3 steps. First frequency is calculated (how much times was the amino acid or nucleotide on that location in the motif), from that you can calculate probability (in the wiki example there were 10 sequences, so each freq. is divided by 10).
Then the log likelihood is calculated, which are the PSSM values. These are in your first matrix (rounded down). The second matrix shows how much the values are relative from your pseudocounts (I assume it was set to default = 0).
Lambda and kappa are estimated to calculate the normalized score (S') for HSP, if you never heard of this before I suggest you read the original psi-blast paper first which you can find here.


text:
Mapping a list of cells in seurat featureplot
I have 209 cells, I clustered them by Seurat to 4 clusters. By Featureplot I am able to track a gene in clusters:

Higher color shows higher expression.
Now, for some genes I want to highlight some cells in Featureplot so that apart from yellow or red colours I want to colour a subsets of cells with another color. I mean I want to map a list of cells in Featureplot or tSNE plot. Let's say I want to know the location of cells 1, 4, 80 and highlight them with another color.
My Seurat object in this link.
Seurat itself beautifully maps the cells in Featureplot for defined genes with a gradient of colours showing the level of expression. Saying I have genes A and B, in excel.
I have coloured cells that express a gene > mean + se, < mean - se or between these values. For instance, for this gene, 36 cells express this gene > mean + se, I want to map these cells in Featureplot or tSNE plot in distinct colour so I can locate them in clusters easily. Something like binary (on off) expression to relative expression.
How I can change the colour of some cells based on my threshold for this gene please?

Do you mean in this exact plot to put a separated color for the cells you want? I don't think this is possible. From what I know, you can put a color code based on gene expression (the situation in your plot) or a color code for a group of cells you want to highlight. But as far as I know both options are incompatible. One thing you can do its to color cells based on gene expression and then highlight the cells you want in a different way. For instance, changing its shape from a dot to a triangle. Do you think this could fit your needs?

Thanks a lot, I think yes. I don't know how to ask Seurat to change the max and min of expression of a certain gene based on my desire. That would be nice if I could color my 209 cells for gene A, if gene A is expression > mean + se, < mean + se, or whatever criteria. The problem is this, in current Featureplot I don't know if a red cell is expressing genes A and B equally, higher or lower.

I will put an answer of how to do it because it is a little bit tricky

To color the TSNEPlot, you can generate a new column in metadata with the expression levels (High, low, etc). Then use pt.shape to set a shape for each identity.
To show binary expression based on expression you first have to define the list of cells that are below or over your threshold. Once you have those lists you can use SetIdent() in Seurat to color those groups.
load(file= "~/seuset_16.RData")

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
seuset_h16 <- StashIdent(seuset_h16, save.name = "expr")
seuset_h16 <- SetAllIdent(seuset_h16, id="res.1")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, pt.shape ="expr")
t

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, colors.use = c("red", "grey"), do.hover =T, data.hover =c("ident", "res.1))
t

Excuse me, I have 4 clusters made by these 209 cells so I want to know the exact location  of cells that are below or over the threshold. For example cell 183 itself places in cluster 0. My question is, by your code is the location of my cells will change in tSNE plot or 4 clusters be reserved???? For instance, this is a link of a picture of my original tSNE plot, I want to know the location of a cell in clusters as it has already been                     https://drive.google.com/file/d/1ipXg2Nb9lM_JgIH6WgVV-PoLhEtzLWif/view?usp=sharing

Sorry, by doing hover I noticed all of cells have been labelled by 1, and high or low ident. That means I can't differentiate the cluster identities anymore (clusters 0, 1, 2, 3). How can color the cells and knowing their cluster identities?

The cluster identity doesn't change, you can add data.hover = c("ident", "res.1"). res.1 is the column name for your clustering at resolution 1 that can be found in object@meta.data. By including this in data.hover you will be able to see the cluster number.

Just edited based on your questions. please upvote answer if this helped. I am new and starting to build reputation! best of luck with your data

Yes and no. Using 3 color schemes won't work but you can leverage the ggplot shapes for this. shape 21-25 allow you to have both fill and color.  do.return = T, returns a ggplot object and using geom_point() allows you to set a shape. You can then use different schemes. for example high threshold open and low threshold filled shapes. What is the scientific information/answer you want to display. There are probably much simpler ways to convey that information.

Sorry, maybe nonsense, I have two lists of genes each have been established to be markers of a developmental stage. you think is there a way to relate my four clusters to these developmental stages based on these two lists?

Since I don't think it is possible to have two color schemes in a single plot (one for expression and one for highlighting some cells) I will suggest other approach.
A possible option to do this could be to change the shape in the plot of the cells that you want to highlight and maintain the color code for expression.
I add a code that works for this purpose with the R object provided. Maybe there are easier ways to do that but this is working:
# Run PCA
seuset_h16 <- RunPCA(seuset_h16,
                     pc.genes = seuset_h16@var.genes)
# Run tSNE
seuset_h16 <- RunTSNE(object = seuset_h16,
                      dims.use = 1:10,
                      do.fast = TRUE)
# Find clusters of cells
seuset_h16 <- FindClusters(object = seuset_h16,
                           reduction.type = "pca",
                           dims.use = 1:10,
                           resolution = 0.6,
                           print.output = 0,
                           save.SNN = TRUE)
# Store output from FeaturePlot (do.return = T)
fp <- FeaturePlot(object = seuset_h16,
                  features.plot = c("DDB_G0273311"),
                  cols.use = c("grey", "orange"),
                  reduction.use = "tsne",
                  do.return = T)
# Add extra column in the data with the cells to be highligthed
fp$DDB_G0273311$data$cellToDisplay <- as.factor(ifelse(row.names(fp$DDB_G0273311$data) %in% c("s1.1", "s1.4", "s1.80"),
"group2",
"group1"))
# Modify color and shape in the plot
# lapply is needed since the object returned in fp is a list
lapply(fp, `+`, geom_point(aes(shape=cellToDisplay, colour = gene, size = 0.4)))

This is the output, with three cells highlighted as triangles and the rest as dots. You can play a little bit with shapes and sizes if you want a clearer display:

First of all thanks a lot for your time, in last piece of code says so                                               > lapply(fp, +, geom_point(aes(shape=cellToDisplay, colour = gene, size = 0.4)))
Error in FUN(X[[i]], ...) : non-numeric argument to binary operator

Is working for me, but try with this: lapply(fp, function(x) { x + geom_point(aes(shape=cellToDisplay, colour = gene, size = 0.4))})

Says that      > apply(fp, function(x) { x + geom_point(aes(shape=cellToDisplay, colour = gene, size = 0.4))})
Error in match.fun(FUN) : argument "FUN" is missing, with no default

I think you miss to copy the l from lapply. We are using the lapply function not apply function. Could be that?

This updates the answer provided by Mack123456 to Seurat 3 and works in 3.1.2. There are some minor difference with colouring, but it is otherwise similar. The HoverLocator plot throws a warning, all about the plot colours and title, but the plot is fine and labels appear as expected.
load(file= "~/seuset_16.RData")

seuset_h16 <- UpdateSeuratObject(seuset_h16)

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(rownames(seuset_h16@meta.data), high)
Idents(seuset_h16, cells = high) <- "high"
Idents(seuset_h16, cells = low) <- "low"
seuset_h16$expr <- Idents(seuset_h16)
Idents(object = seuset_h16) <- "res.1"

t <- DimPlot(seuset_h16, pt.size = 2, shape.by = "expr")
t

t2 <- DimPlot(seuset_h16, pt.size = 2, cols = c("red", "orange","green","blue"))
t2 <- HoverLocator(plot = t2, information = FetchData(object = seuset_h16, vars = c("ident", "res.1")))
t2


text:
Running python script from terminal: umap issue
Some time ago I asked the question:
RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'
It was solved by following the answer in it, however, it is still not working when I run a python script from terminal, giving the same error:
python python/umap.py -path
"/Users/nikitavlasenko/python_scripts/output/" -data
"/Users/nikitavlasenko/python_scripts/data_files/WT_SCID_CCA.csv"

Traceback (most recent call last):
    File "python/umap.py", line 26, in
      result = umap.UMAP().fit_transform(data.T)
  AttributeError: module 'umap' has no attribute 'UMAP'

From inside spyder the same script runs perfectly, but somehow when running it from the terminal it is not successful. I am new to python, so maybe I am not pointing correctly to umap library? Maybe it should be imported in a different way? Here is the whole script:
import umap
import numpy as np

def getopts(argv):
opts = {}  # Empty dictionary to store key-value pairs.
while argv:  # While there are arguments left to parse...
    if argv[0][0] == '-':  # Found a "-name value" pair.
        opts[argv[0]] = argv[1]  # Add key and value to the dictionary.
    argv = argv[1:]
return opts

if __name__ == '__main__':
  from sys import argv
  myargs = getopts(argv)
  print(myargs)
  path_to_save = myargs['-path']
  path_to_data = myargs['-data']
  data = np.loadtxt(open(path_to_data, "rb"), delimiter=",", skiprows=1)
  result = umap.UMAP().fit_transform(data.T)
  np.savetxt(path_to_save + "/umap.csv", result, delimiter=",")

Python version that is used by spyder and in the terminal seems to be the same: 3.6. I am using Mac OSX.

Add a print(umap.__file__) early on. It's likely that you're getting a different copy of umap.

You answered your own question, yet it is not a coincidence that you have to change the name of your script.
When Python executes your custom script, it reads:

import umap

Python will first search the local sources (the directory of the script) , and afterwards the installed sources.
In this case Python identifies the local source to be umap.py, and imports the (currently executed) script, while foregoing the installed sources.
Since you are calling functions, not defined in your current script, Python will not be able to fully execute your script.
For more information on importing see Module search path at docs.pyhton.org

It is pretty weird but the only thing that I needed to do is to rename my python script from umap.py to umap_gen.py and right away it started to work. So, I guess there can be such name collisions which we should avoid. Although I would prefer the file to be named as umap.py still...


text:
Problem installing last build of Sniffles
I have followed the steps described here to install Sniffles on my Mac Book Pro running macOS High Sierra.
However, I am struggling with this step:
cmake -D CMAKE_C_COMPILER=/opt/local/bin/gcc-mp-4.7 -D CMAKE_CXX_COMPILER=/opt/local/bin/g++-mp-4.7 ..

I do not have an /opt folder, so I get the following output:
-- The C compiler identification is unknown
-- The CXX compiler identification is unknown
CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_C_COMPILER:

    /opt/local/bin/gcc-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CC" or the CMake cache entry CMAKE_C_COMPILER to the full path to
  the compiler, or to the compiler name if it is in the PATH.

CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_CXX_COMPILER:

    /opt/local/bin/g++-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CXX" or the CMake cache entry CMAKE_CXX_COMPILER to the full path
  to the compiler, or to the compiler name if it is in the PATH.

-- Configuring incomplete, errors occurred!
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeOutput.log".
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeError.log".

Unsurprisingly, the error file says that /opt/local/bin/gcc-mp-4.7 does not exist.
So I tried:
cmake -D CMAKE_C_COMPILER=/usr/bin/gcc -D CMAKE_CXX_COMPILER=/usr/bin/gcc ..

I get another long message but the most important bit is:
CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

I do not know what to do to solve this error.
I do have gcc here: /usr/bin/gcc and cmake version 3.8.1.
I have created an issue on the Sniffles GitHub repository.
Update 13th June 2018
I installed the latest version of gcc:
brew install gcc
Warning: gcc 8.1.0 is already installed and up-to-date
To reinstall 8.1.0, run `brew reinstall gcc`

Then, I tried to make cmake use gcc and g++ in the /usr/local/bin folder.
cmake -D CMAKE_C_COMPILER=/usr/local/bin/gcc-8 -D CMAKE_CXX_COMPILER=/usr/local/bin/g++-8 ..

This recreated the same problem:
CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

Do I need to set the OpenMP_C*_FLAGS to something?

Do you have OpenMP installed? It's required.

@DevonRyan OpenMP is installed as part of GCC. The problem is that OP isn’t using GCC.

You need to install GCC, or a newer version of Clang.
Your compiler — though the binary is called gcc — is actually Clang, and the Clang version that’s shipping with macOS/Xcode still does not support OpenMP.
The easiest course of action is probably to install GCC via Homebrew (brew install gcc) and then configure Sniffles to use Homebrew’s GCC (i.e. setting the CMAKE_C*_COMPILER paths to /usr/local/bin/g{cc,++}-$version).

I have tried what you said, unsuccessfully. I have updated my original question.


text:
Can you query the RCSB PDB API with a ligand and return the IDs of all entries in which it appears as a free ligand?
If you're looking at the entry for a ligand, for example, adenosine, there's a link on the webpage to entries for which this entry exists as a free ligand.
However, I want to access this via the API and can't seem to find docs/info on how to do that.
I have been looking at the fields you can use to get custom reports.

I'm not sure about the RCSB, but you can get such a list from the PDBe REST API, as documented in the compounds section:
http://www.ebi.ac.uk/pdbe/api/pdb/compound/in_pdb/ADN
Alternatively, from an older service:
http://www.ebi.ac.uk/pdbe-srv/pdbechem/PDBEntry/download/ADN


text:
determine if ChIP-seq peaks are broad or narrow
Is there a method to determine if the peaks are broad or narrow? ENCODE provides some guidelines:

Although those cover common histone marks, there are many others. If you are using one of the ones that are not listed, is there a way to determine the appropriate peak width?
This question was inspired by this previous post.

I'd say mostly this is a question of understanding the underlying biology and the relevant literature.
If it is not known in the literature whether a mark is peaky or broad, evidence might come from FISH studies or Low throughput qPCR. Another way to look would be to examine the signal expressed as fold enrichment over input on a genome browser and look if peaks tended to be small (a few hundred bp, e.g. a transcription factor), or large (several 1kb, e.g. H3K27me3) or contain both (e.g. RNA polymerase II).

That makes sense, but H3K27ac is frequently thought to be broad and the peaks look broad, but is apparently narrow.

Always these things are going to be a judgement call. There are going to be marks that some people count as braod that others count as narrow, and vice versa.


text:
How can I determine a mean sequence divergence for 10k sequences?
I am trying to analyze a number of repetitive sequences and as one step want to calculate a sequence divergence between the elements I found.
Now in theory I wanted to generate a MSA of the sequences (~100bp times ~15k) and use the MEGA7 option to calculate pairwise and mean divergence. But sadly I run into the problem:

Error# 4511 encountered: No common sites were found for the sequence
  pair (193,4).

and my MSA I generated with MAFFT looks like:
>h1
-------gat-tt------------tgta--t---c-t------a---ta---t--t-a-
g-ca--tt-c-tt-t--at----------c----------a--------------ttt--
--------a---t-aata----------c---------c-------------------t-
------ag-a--a-a--------------------t-------------g----------
----t------tt-t-a-a---aa------------------------------------
------t--g--------c-----a-a-t-----t-----t-------------------
-------------a-------a--a--c---c-----a---t-t------------tt--
------------------t---g---g---------t------------------t-t--
---t-----g-a--tc--------------------at-------a--gt-ag----gc-
---t------a---t--t------t-t---c--a--a----a---a--aaatatc---

basically my sequence normally is 95% conserved just single individuals lead to the accumulation of those gaps, I am not sure how I can reduce these.
How can I calculate the sequence divergence? Any other programs? Any other methods?

consider removing those sequences? It's unlikely that any numbers you obtain for them will be very meaningful.


text:
CNSeg no method coercing S4 class to vector
I'm attempting to analyze some genomic data with the help of a package called CNSeg. The package is located here and the documentation is here.
My input data looks like:
> str(segments)
'data.frame':   11897 obs. of  7 variables:
 $ X                : int  0 1 2 3 4 5 6 7 8 9 ...
 $ SegmentID        : int  72 73 74 75 76 77 78 79 80 81 ...
 $ Chromosome       : int  1 1 2 2 2 3 3 3 3 3 ...
 $ StartPosition    : int  754192 145260908 21494 141215321 141275624 63411 69812903 69884262 126473310 126790130 ...
 $ StopPosition     : int  145258178 249212878 141214996 141275051 243052331 69811900 69884106 126457276 126772699 197852564 ...
 $ Median.Log2.Ratio: num  -0.014 0.311 -0.003 0.059 -0.012 -0.018 -0.106 0.007 -0.171 0.001 ...
 $ FileName         : Factor w/ 95 levels "TSB02183","TSB02184",..: 1 1 1 1 1 1 1 1 1 1 ..

(Note: I did import the segments data with FileName as a vector instead of a factor as well. No difference in the error below).
I followed the vignette to the letter and, at the end of the vignette, the creator tells me that I can use the dist() method. When I try to use it, I get an error: Error in as.vector(data) :
  no method for coercing this S4 class to a vector.
Here is my session information:
> segments <- read.csv("Probe_Segments_CN.csv")
> cnseg <- CNSeg(segList = segments, chromosome = "Chromosome", end = "StopPosition", start = "StartPosition", segMean = "Median.Log2.Ratio", id = "FileName")
> rdseg <- getRS(cnseg, by = "region", imput = FALSE, XY = FALSE, what = "mean")
Processing samples ... Done
> data("geneInfo")
> geneInfo <- geneInfo[sample(1:nrow(geneInfo), 2000), ]
>
> rdByGene <- getRS(cnseg, by = "gene", imput = FALSE, XY = FALSE, geneMap = geneInfo, what = "median")
>
> reducedseg <- rs(rdseg)
> f1 <- kOverA(5, 1)
>
> ffun <- filterfun(f1)
>
> filteredrs <- genefilter(rdseg, ffun)
> filteredrs <- madFilter(rdseg, 0.8)
> dist(filteredrs)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] tools     stats     graphics  grDevices
[5] utils     datasets  methods   base

other attached packages:
[1] CNTools_1.34.0    genefilter_1.60.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17         AnnotationDbi_1.40.0
 [3] BiocGenerics_0.24.0  splines_3.4.0
 [5] IRanges_2.12.0       bit_1.1-14
 [7] lattice_0.20-35      xtable_1.8-2
 [9] blob_1.1.1           parallel_3.4.0
[11] grid_3.4.0           Biobase_2.38.0
[13] DBI_1.0.0            survival_2.41-3
[15] bit64_0.9-7          digest_0.6.15
[17] Matrix_1.2-9         S4Vectors_0.16.0
[19] bitops_1.0-6         RCurl_1.95-4.10
[21] memoise_1.1.0        RSQLite_2.1.1
[23] compiler_3.4.0       stats4_3.4.0
[25] XML_3.98-1.11        annotate_1.56.2

I was not able to get the CNSeg clustering method to work, but the work-around was to export the data and cluster it independently of the CNTools package (thanks to Llopis's answer below). I'll post the analysis below.

require(CNTools)
segData <- read.csv("result_cnv.csv", stringsAsFactors = FALSE)
head(segData)

# Create inital CN object
cnseg <- CNSeg(segList = segData, chromosome = "Chromosome", end = "StopPosition", start = "StartPosition", segMean = "Median.Log2.Ratio", id = "FileName")
cnseg

# Create inital RD object
rdseg <- getRS(cnseg, by = "region", imput = FALSE, XY = FALSE, what = "mean")
rdseg

# Collect gene information
data("geneInfo")
geneInfo <- geneInfo[sample(1:nrow(geneInfo), 2000), ]

# Create an RD based on gene information
rdByGene <- getRS(cnseg, by = "gene", imput = FALSE, XY = FALSE, geneMap = geneInfo, what = "median")

# Initalize reduced segment
reducedseg <- rs(rdseg)

# Create a function that evaluates to TRUE if at least 5 of the argument elements are larger than 1
f1 <- kOverA(5, 1)

# Create a filter based on f1
ffun <- filterfun(f1)

# Use the CNTools genefilter
filteredrs <- genefilter(rdseg, ffun)

# Use the CNTools madFilter
filteredrs <- madFilter(rdseg, 0.8)

filteredrs

# Write the filtered data to file
write.csv(attributes(filteredrs)[1], "CNseg_CNOut.csv")

# Reading the recently-written file helps with formatting
CNseg <- read.csv("CNseg_CNOut.csv")

# Keep only the columns with actual sample values
CNseg <- CNseg[5:ncol(CNseg)]

# Calculate eluclidian distance between samples
d <- dist(t(CNseg), method = "euclidean")

# Calculate clustering
hc1 <- hclust(d, method = "ward.D")
plot(hc1, cex = 0.6, hang = -1, main = "Clusters of Copy Number Alterations", xlab = "Euclidean Distance")

If you are using the latest version of Bioconductor you should always look at the latest documentation of the package in Bioconductor. There could be changes from 2011 to 2018 (In this case there aren't any changes in the section you mention).
The problem is the package doesn't define how dist should behave with that object class (filteredrs is a RS object if I read the source code well).
You can extract the metrics to calculate the distance yourself or ask the maintainer to add this method to the class. However it could be that the maintainer answer you if you post the question in the support site of Bioconductor.

Good catch--I'll look at more updated documentation in the future and I'll ask the question on the Bioconductor support site. And yes, you read the code correctly I believe. How would you extract the metrics? It's in the values pane as a "Large RS" value with a dataframe inside, I just don't know how to get to the dataframe.

EDIT: Solution is attributes(filteredrs)[1], since "$" doesn't work on the object.

@CalendarJ You can post (and accept) your own solution, that way you can help future readers (comments are harder to find [not indexed by google] and can be deleted).


text:
Public access to genomics databases
I'm a statistician, and am interested in applying the theory of algebraic Markov models to genomics. Here's one paper I'm interested in: Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. I'm also interested in other diseases (i.e. not just bipolar).
Are there any genomics databases that I can access publicly?
There are some links in Human Connectome Project, but I'm not sure if these are genomics databases (or which papers have used their research).

You're going to struggle to find public databases with genotype data for individuals because of concerns about de-identification. You can request access to the server where the data from that paper is stored here. If it's possible to run your analyses on summary statistics, these can usually be downloaded without restrictions.

Could you clarify what kind of data you need? "genomics databases" is a very broad term and most of what I would call genomics data is publicly available. It would be helpful if you could [edit] your question and explain exactly what sort of dataset you would want to work on.

Most genomics data that is published should be available in GEO (or SRA). Usually in the paper they mention the accession number. However, not every study ends up in there, and sometime studies that are in GEO are not reproducible (poor description of the samples, or design). But if you look at the original papers, most of the times you'll find an accession number.

Another places you could look for SNP data would be 1000 genomes, the ExACT database and the genomAD database. Cancer specific data is available at the GDC and ICGC data portal.
Unfortunately a lot of what you'd like to work with is probably covered by restricted access due to privacy concerns. Lots of this data is stored in the dbGap database and the EGA. If you are an academic or employed by a company with a legitimate interest in the data you can usually apply for access to the restricted data, and I've never known anyone to be turned down if you have a legitimate reaons to access the data and can demonstrate your computer systems meet the required security standards.

Yet another place I find useful is provided by Repositive.
You can think of it as a meta-search engine which includes most of the databases already stated and links hits to relevant papers.
With respect to your question about bipolar, this search could get you started.


text:
Violinplot of gene expression
I have plotted the log normalized expression of two genes by violonplot for 4 clusters. I have links to my pictures and Seurat object too. My problem is this; in violin plot I can not see the mean or any centennial tendencies so that I don't know if two genes is expressing higher or lower in contrast to each other in each cluster.
violonplot vlnPlot(seuset,c(""))
Actually when I am thinking I see I need a violin or box plot showing  the  ratio of these genes to each other in each cluster and how nice to have mean also in the plot
Somebody knows how I can do that?

Please, remember to add the code you use to make it easier to provide the accurate advise to help you. But do you want to see the mean of the cluster or to see the differences of genes between clusters?

You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.
(1) First, notice that vlnPlot() is deprecated. Use VlnPlot().
You can try using the parameter do.sort=T:
VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)
Alternatively, you can return the ggplot2 object and then plot the means etc:
p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p

(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.
In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.
Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group.
(a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.
Here are some functions for retrieving and plotting data from the object:
ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=1)

stripchart(datatoplot2,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=2)

excuse me, with this command i have this picture in the link for my four clusters p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)   https://drive.google.com/file/d/1r6eGQB225_jwtf7AWQo6Z2FKj4eVro42/view?usp=sharing, how i can have the plot for overall data side by side this plot?

i plotted that for all of cells but i don't know how to make a 5 violin together. https://drive.google.com/file/d/1MarsjXbTf0jg8e8e-1MTARNFQsOC9svw/view?usp=sharing   so i would be able to compare each cluster mean with the overall  mean, how i can have the plot for overall data side by side this plot? i plotted that for all of cells but i don't know how to make a 5 violin together


text:
Spearman correlation for large dataset
I have two datasets (DataA and DataB) and I want to find the Spearman correlation between genes and also pull out the gene names (stored in first column of dataset) in R. I am using fread from read.table to read the file and cor.test to find Rho and p-value.
The code works fine for a smaller data set of 10 obs, 10 variables. But, my actual data set is larger (DataA - 795 obj, 542 variables and DataB - 925 obs, 542 variables).
When I try to run the code on my actual data, it is consuming a large amount of time (it has been 23 hrs since my code is running). Is there a way to optimize this code? Or is there a bug in the code?
library(data.table)
A.data <- fread("DataA.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)
B.data <- fread("DataB.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)

correlation.result <- data.frame()

## Computing Correlation

for (i in 1:nrow(A.data)) {
  for (j in 1:nrow(B.data)) {
       correln <- cor.test(as.numeric(A.data[i,2:ncol(A.data)]), as.numeric(B.data[j,2:ncol(B.data)]), method = "spearman", exact = FALSE)
       rho <- correln$estimate
       p.val <- correln$p.value
       Gene <- paste(A.data[i,1], B.data[j,1], sep="_", collapse=" ") # To pullout gene name pairwise
       comb.data <- data.frame(rho, p.val, Gene)
       correlation.result <- rbind(correlation.result,comb.data)
 }
}

Your code looks fine, but rbind to progressively build and resize a data frame can be gluttonous in terms of memory and time, especially given the A-row x B-row times you are doing this operation, which is fine for very small datasets, but not often realistic for real-world data.
Instead, perhaps cat results to a text file, where results are written in tab-delimited columns. You can then read in that text file later on.
In addition to serializing or saving your results to text, which make them easier to work with in non-R contexts, you won't have to deal with R's memory allocation problems when running rbind repeatedly. You just write out the answer to a file and read it in later on, in one pass.
Another approach is to preallocate memory to a presized data frame, and then write to that allocated space.
You know the size of your inputs ahead of time, which makes this approach tenable. In this case, it looks like you are storing two floats and a string (rho, p.val, and Gene) and that you are going to store A-row x B-row number of them.
Here's one way you might go about this:
# create a presized data frame
Arows <- nrow(A.data)
Brows <- nrow(B.data)
correlation.result <- data.frame(t(matrix(c(0.0, 0.0, "foo"), nrow=3, ncol=Arows * Brows)), stringsAsFactors=FALSE)
colnames(correlation.result) <- c("rho", "p.val", "Gene")

# assign values to that presized data frame
Acols <- ncol(A.data)
Bcols <- ncol(B.data)
for (i in 1:Arows) {
  # get A data
  Ai <- as.numeric(A.data[i, 2:Acols])
  Aname <- A.data[i, 1]
  for (j in 1:Brows) {
    # calculate row offset in presized data frame
    ij <- (i - 1) * Brows + j
    # get B data
    Bj <- as.numeric(B.data[j, 2:Bcols])
    Bname <- B.data[j, 1]
    # run test
    correln <- cor.test(Ai, Bj, method = "spearman", exact = FALSE)
    rho <- correln$estimate
    p.val <- correln$p.value
    Gene <- paste(Aname, Bname, sep="_", collapse=" ")
    # assign value to row of result data frame
    correlation.result[ij,] <- c(rho, p.val, Gene)
  }
}

This avoids R needing to resize and recopy the data frame contents as they grow, which I suspect is the major reason it is taking so long to run.
Chapter or "Circle" 2 of Patrick Burns' book The R Inferno discusses the problems with growing objects with rbind and ways to handle it — you might find this a useful read. This book is a PDF download via your search engine of choice.

Is there a reason for the downvote(s)?

I think it is bad for a community when downvoting without any explanation. Please explain yourself.

You can use corr.test() from psych package, e.g.
ct <- corr.test(t(A.data), t(B.data), method='spearman', adjust='none')

then ct$p and ct$r would give you (uncorrected) p-value and correlation matrices of size 795 x 925. Note that t() is required to compare matrices row-wise.
That would be orders of magnitude faster than loop-based approach.


text:
Can I ask DESeq2 which variable alone explains which gene's behavior the best?
Imagine I have a dataset of highly periodic gene expression, taken at densely spaced timepoints. I know that the expression is periodic, but different genes have different offsets.
It is trivial to convert my timestamps into a cosine (cos_time = cos(offset + time*pi*2/PERIOD)), so that this new variable is linear with the genes' expressions -- but only with the ones that have a peak at offset.
What if I could have multiple cos_times (let's call them phases) and then ask DESeq2: "OK, which gene's behavior is most linear with which phase?" Barring comparing p-values of each LRT test, is there a better way?
Disclaimer: I know about cosinor, Lomb-Scargle periodograms, ARS, JTK etc, and I know they are likely more applicable in this particular case. I'm interested in whether it is in general possible to ask such a question from DESeq2, even if the variables in question are not periodic. Phase just seemed like the simplest concept to demonstrate this on.

If you have densely spaced time points the using non-linear regression to directly fit the signal to a cosine or sine curve, which would directly produce the offset, would be much simpler. DESeq2 wouldn't be the appropriate tool then, you'd use base R and the nls() function.

If I were to do this, my first guess would be to add in the phases as a separate factor in the statistical model, something like:
~ Condition + sex + phase1 + phase2 + phase3 + ...

And then run LRTs removing phases from the model one at a time? This will leave me with p-values from individual tests, and comparing p-values to each other is probably going to put the wrath of statisticians upon me...
Or is there a way to run DESeq2 in such a way that it spits out the "best phase" for each row?


text:
Better aligner than bowtie2?
Bowtie2 is probably the most widely used aligner because of it's speed. Burrow-wheeler (BW) algorithms (including bwa) tend to be faster. However, they have limitations when it comes to aligning very short reads (e.g. gRNA). Also, setting maximum number of mismatches allowed is complicated by the seed length, overlaps and other parameters.
I wonder if there is any better multi-purpose aligner out there. May be with algorithm other that BW. One which allows special cases e.g. allowing shord reads and high number of mismatches.

How long is "very short"? Bowtie2 is not optimized for <=36bp query sequences; nor most aligners developed in the past 5 years or so. Bowtie1 or bwa-aln will work better.

@user172818 When you say "not optimized", do you mean using the default settings, or more generally ?

@bli in general. Aligning 100bp reads is a very different problem from aligning <36bp reads.

If I have learned one thing working with bioinformatics over the years, that is there is no such thing as a "multi-purpose" or "overall best" tool out there. Different algorithms are good for different types of data. My recommendation would be to arrange your workflow in a modular way and plug in the best aligner for the type of data you are working on, rather than relying on a "best" solution.

As of 2020-09-09, according to google scholar, the 2012 bowtie2 paper has 21647 citations, the 2009 bwa paper has 25045 citations. It may be true that bowtie2 is the most widely used aligner, but not probable.

Bowtie2 is no longer the fastest aligner. Salmon and Kallisto are much faster, but have been designed to optimise RNASeq mapping. Their speed is gained from avoiding a strict base-to-base alignment, but they can output mostly-aligned reads (i.e. position-only, without local alignment) as pseudo-alignments. See here for more details.
Both Kallisto and Salmon can do additional bootstrapping that is interpreted by sleuth for improved performance in isoform detection. They can also output counts that are equivalent to read-level counts from other programs, which can then be used by other downstream gene-based differential expression analysis software (e.g. DESeq2). Salmon has additional options that can  correct mappings for sequence-level and GC bias.
HISAT2 is from the same group as Bowtie2, and does the same sort of stuff, but with a few optimisations added on top. In particular, it's much better at working out split reads from RNASeq runs. Like Bowtie2, it will do local alignment of reads.
For quick alignment of long reads, minimap2 works well. For high-accuracy alignment (but comparatively slower), LAST works well.
Most bioinformaticians seem to prefer STAR for things that Bowtie2 was previously used for. I'm not yet convinced it's a better alternative, and currently prefer HISAT2 for high accuracy short-read alignment.

Seconded for kallisto, it is very good. But currently downstream support is limited to the other tool from the pachter lab for differential gene expression analysis, called sleuth.

One advantage to STAR aligner in Bioinformatics and the reason it is used at my work, is the possibility of reading Chimeric alignments (for detecting e.g. circular RNA through custom coding).

Salmon and Kallisto aren't really aligners and are only useful for RNA-seq, though they are great and very quick. @NatWH Kallisto can also be used with DESeq2 or edgeR.

@Jared Andrews at least when I last used it, Kallisto only outputs TPM directly, whereas the normalisation applied by DESeq2 requires unnormalised counts. Are there settings to change the output of Kallisto?

@NatWH It outputs estimated counts as well. And you can see it mentioned as usable for DESeq2 in the DESeq2 manual.
From kallisto manual: "abundances.tsv is a plaintext file of the abundance estimates. It does not contains bootstrap estimates. Please use the --plaintext mode to output plaintext abundance estimates. Alternatively, kallisto h5dump can be used to output an HDF5 file to plaintext. The first line contains a header for each column, including estimated counts, TPM, effective length."

@Jared Andrews, yes it is possible by summing over transcript estimated abundances to produce integer-like counts, but Lior Pachter argues here that this should not be done. What are your thoughts?

@NatWH Basically the same as the response to his comment. Important to note that those concerns were raised nearly 3 years ago and both salmon and kallisto have been updated several times since then. Regardless, since it's recommended by the DESeq2 authors themselves, I am inclined to trust them - see their paper on the subject. And you don't have to summarize to gene-level counts if you don't want to/don't trust it.

@Jared Andrews alright. I'll leave my comments up in case any future visitors benefit from the links which have been posted.

I second what @JaredAndrews said: Salmon and Kallisto are not aligners or direct alternatives to Bowtie. I can't how one would replace the pipeline Bowtie mapping -> macs + bigwig with those, or pretty much anything related with sRNA analysis. A read aligner will output the mapping position of a read, whereas salmon and kallisto will output summarized estimates of gene expression based on the read + transcript properties. Since the OP does not state gene / isoform expression as a use case, I have down-voted this answer because of that paragraph.

Kallisto can output position. What it doesn't do is a careful local alignment (and the question didn't specify whether local alignment was important). I indicated other alternatives where local alignment is important.


text:
How do I estimate branch length of a maximum parsimony phylogeny in R?
I am trying to estimate a phylogenetic tree in R with ape and phangorn. With pratchet() and optim.parsimony(), I can estimate a phylogenetic tree, but branch lengths are not estimated. Is there a way to estimate a phylogeny with branch lengths with the following made-up data:
x <- data.frame("L1" = c(0,0,1,0,0,0,1,1,1,1,0,0,0,0,1), "L2" = c(0,1,1,0,1,1,0,1,1,0,1,1,0,1,1),"L3" = c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1), "L4"=c(0,0,1,1,0,0,0,0,1,1,0,0,0,0,1))

Don’t quote your identifiers! It’s L1, L2, etc.; not "L1", "L2" …. R unfortunately allows is but this is a language flaw that should not be exploited.

Can you provide a link to further information @konrad-rudolph

@NathanS.Watson-Haigh Unfortunately I’m not sure what further information you’re looking for, I’m afraid.

@Nathan S. Watson-Haigh's answer does give branch lengths, but NOT parsimony branch lengths, which may be what the OP wanted.  nnls.phylo is giving distances calculated from a pairwise distance matrix. To get the parsimony branch lengths (i.e. number of changes mapped on to each branch), score the tree against the dataset using acctran:
treeRatchet <- pratchet(primates)
treeRatchet.BL <- acctran(treeRatchet,primates)
plot(treeRatchet.BL,main='Phylogeny (Parsimony branch lengths)')

Note that the branch lengths in the parsimony phylogram are 'chunkier' than those in the nnls.phylo phylogram -- that's because

they are all integers, and
they aren't corrected with reference to a model of substitution that considers multiple hits or base composition.

Terminology
First, a couple of corrections to your use of terminology. A tree must have branch lengths associated with it to be called a phylogenetic tree as it represents both the relationships between OTU/taxa/leaves AND the distances between them (i.e. the branch lengths. A tree where branch lengths have no meaning is called a cladograms or a dendrogram. You can tell if a tree is a cladogram since all the OUTs/taxa/leaves end at the same position and internal branch lengths are usually (not always) of equal length. Whereas a phylogeny the OTUs/taxa/leaves will typically look "staggered" and internal branches will be of different lengths.
Therefore, If the tree you have estimated, truly doesn't have branch lengths estimated and associated with it, it is actually a cladogram.
Example Derived from phangorn Vignette
library(phangorn)

# Get example data
fdir <- system.file("extdata/trees", package = "phangorn")
primates <- read.phyDat(file.path(fdir, "primates.dna"), format = "phylip")

# Calculate pairwise distances between OTUs using Maximum Likelihood
dm <- dist.ml(primates)

# Use the distances to cluster the OTUs using UPGMA algorithm
# This creates a dendrogram
treeUPGMA <- upgma(dm)
plot(treeUPGMA, main="UPGMA")
dev.off()

# Search for different, similar trees to the initial UPGMA tree
# looking for a tree with the lowest parsimony score
# There are 2 search algorithms implemented:
#   optim.parsimony()
#   pratchet()
# These are still cladograms
treePars <- optim.parsimony(treeUPGMA, primates)
treeRatchet <- pratchet(primates, trace = 0)
plot(treePars, main="optim.parsimony()")
plot(treeRatchet, main="pratchet()")
dev.off()

# Calculate branch lengths using the most parsimonious tree identified by
# the pratchet() searchh algorithm
phylogeny <- nnls.phylo(treeRatchet, dm)
plot(phylogeny, main="Phylogeny")
dev.off()

Choices, Choices, Choices
As you can see from the above example, there are lots of steps to calculating a phylogenetic tree. At each of these steps there may be multiple algorithms to choose from. You choice will depend on the type of data you are analysing, the size of the data set, your available compute resources and the algorithms implemented and available through the tools you use.
Ultimately, you will have to know a fair bit about these algorithms to make an informed choice.


text:
Is it possible to identify cells that are expressing two or more genes in Seurat?
I'm interested in looking into cells that are positive for two (or in some cases more) genes. I know I have some double positive just by looking at the FeaturePlot of those genes, but now I'm trying to figure it out how many are double positive and further percentage wise. Example:
50 cells are positive for Myh11
100 cells are positive for Acta2
30 cells are double positive
60% of Myh11+ cells are Acta2+
30% of Acta2+ cells are Myh11+
Is the data used to generated the FeaturePlot the scale.data? If so, I could just make a new dataset out of that, but I'm not sure if that is the case.

The counts stored in the Seurat object are: raw counts (seuratobject@raw.data), the log + normalized counts (seuratobject@data), and the scaled counts (seuratobject@scale.data). FeaturePlot() plots the log + normalized counts.
In order to identify double-positive cells, you need to identify cells that express a gene (i.e. positive for a gene) and that is not trivial: you need a cutoff value for each gene, and also to distinguish between true and dropout zero counts.
Provided you solve these, you can retrieve your cells by:
MYH11.cutoff <- 1
ACTA2.cutoff <- 1
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff))
length(which(FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff & FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))

you can substitute the length(which( by just using sum: TRUE is evaluated as 1 and then you end up with the same result (but a bit faster).

Sorry   what  would be the code if   we have    3 clusters of  cells and we we want to show the cluster in which the ratio of gene named DDB_G0267412 is two times higher than DDB_G0277853?


text:
Download proteomes from NCBI based only on binomial names
I'm trying to download around 2000 proteomes from NCBI, but all I have is the genus and species e.g: lophodermium_seditiosum trichoderma_harzianum
Does anyone know a good way of doing this in script form?
My code at the moment is:
from Bio import Entrez
from Bio import SeqIO
import os

file_name = "proteome1.faa"
net_handle = Entrez.efetch(db="protein", id="trichoderma_harzianum", rettype="gp", retmode="text")
out_handle = open(file_name, "w")
out_handle.write(net_handle.read())
out_handle.close()
net_handle.close()
print("Saved")

for seq_record in SeqIO.parse("proteome1.faa", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

But I think it wants a number as the id, I don't know if it's possible to convert species name to a numerical ID within the script, or if that's the best way of going about this.

The general procedure is explained here:
First step is to search the Taxonomy database with the organism name. Accepted common names usually work at all taxonomic levels. Use the scientific name or formal name if no results are obtained with the common name.
In code:
from Bio import Entrez

Entrez.email = 'user@example.org' # Put your email here

def get_tax_id(species):
    species = species.replace('_', '+').strip()
    search = Entrez.esearch(term=species, db='taxonomy', retmode='xml')
    record = Entrez.read(search)
    return record['IdList'][0]

if __name__ == '__main__':
    organisms = ['lophodermium_seditiosum', 'trichoderma_harzianum']
    taxids = [get_tax_id(organism) for organism in organisms]
    print(taxids)

This gives us the taxids for our species, respectively 128023 and 5544. Now we can use these to get the proteomes. For Lophodermium seditiosum we find two results, but for Trichoderma harzianum there are 65900 results, so you might want to tweak the query to filter the results. For example, if you only want the sequences appearing in PDB you add AND pdb[filter] so that the query becomes txid5544[Organism:exp] AND pdb[filter].
def get_proteomes(taxid):
    query = 'txid{}[Organism:exp]'.format(taxid)
    handle = Entrez.esearch(db='protein', term=query)
    result = Entrez.read(handle)
    if int(result['Count']) <= 2000:
        ids = ','.join(result['IdList'])
        proteomes = Entrez.efetch(db='protein', id=ids, rettype='fasta').read()
        return proteomes
    else:
        raise RuntimeError('Too many results!')

if __name__ == '__main__':
    print(get_proteomes('128023'))

Returns:
>ADI44294.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

>ADI44293.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

So I'm getting a 'Supplied ID parameter is empty' error from this code, the functions are defined at the beginning and are the same as in the answer
`listoforganisms = [x.split('\t')[0] for x in open("OGTlist.csv").readlines()]
IDlist = []
#convert to id
if name == 'main':
organisms = listoforganisms
for organism in organisms:
taxid = get_tax_id(organism)
IDlist.append(taxid)
print (taxid)
print(get_proteomes((taxid)))
`

Can you post a small sample of your OGTlist.csvfile?

Got it working with, `if name == 'main':
organisms = listoforganisms
for organism in organisms:
taxid = get_tax_id(organism)
print(organism)
print (taxid)
stringid = str(taxid)
    print(get_proteomes(stringid.strip("'[]'")))`

Just needed to strip the list elements to the right form

You can use the NCBI Datasets to download the protein sequences of all annotated protein-coding genes for an organism starting with the organism name. There is a web-interface that allows exploring by taxonomy or specific organism search as well as a command line utility and an API. For example, the command line utility can be used as follows:
$ datasets download genome taxon 'trichoderma harzianum' \
    --annotated --exclude-genomic-cds --exclude-gff3 \
    --exclude-rna --exclude-seq --reference
Collecting 1 genome accessions [================================================] 100% 1/1
Downloading: ncbi_dataset.zip    3.92MB done

This downloads a zip archive containing the protein.faa file for the Trichoderma harzianum reference assembly.
Note, unlike the Entrez method described above, this will not download all proteins for a given organism from the NCBI Protein portal. Thus, Lophodermium seditiosum does not return any results as there is no available genome for this assembly in the NCBI Assembly database.
On the other hand, the NCBI Datasets method will be significantly faster and return cleaner results. I would argue that if you are interested in obtaining the 'proteomes' of a given set of organisms, this is a better choice as it includes only those cases where a complete proteome is available based on the annotation of the genome and does not include cases where only a few proteins from an organism are available.


text:
The slopes of gene networks
I have a gene network that I clustered into 8 parts. I calculated the degree distribution of the clusters using Cytoscape and the slope of the clusters range from $0.7$ to about $1.3$.
From literature (reference) usual scale free networks (real world networks) have a slope range of $2>γ>3$. What could be a reason for my network to have the smaller slope value?
Edit: The number of nodes of each cluster ranges from 431 to 187 and the number of edges range from 1009 to 263. The network is bipartite in nature.  The clusters are all very sparse which can be understood by the number of edges in the largest and smallest cluster. Are real world biological networks not always scale free. Can anyone point me towards publications that mention this fact. I am just confused about why my network doesn't conform to the norm but maybe it doesn't need to.

Could you add more information about your network? How many nodes does it have? How are them distributed in your 8 clusters? What kind of data are you using to create it (It seems like gene expression, but I'm not sure) ?

Yes, please [edit] and explain more about your network. 8 clusters suggests your network is likely very small. Is it also sparse, as biological networks often are? You also seem to be starting with the assumption that your network is scale free. Is it? It might not be, after all.

@terdon Aren't real world biological networks scale free in nature? Or that's what I have been led to believe. I have not been able to find any publications to the contrary yet. What more information do you think I can provide about the network. I have mentioned the number of nodes and edges.

Actual networks are rarely scale free.

@Devon, this seems to be the answer rather a comment. Thanks for providing the reference (I lost it/I didn't save it when I read it). @ The Last Word: Despite not being a scale-free you can still work with your clusters, but take note that the topology of your network is not to be trusted.

@Llopis I'm hoping that someone else wants to give a more thorough answer than my link :)

@Devon I think there was another paper about networks in biology, but in general it seems like people is using topologies analysis assuming too many things without (first) checking them.


text:
Why does `pratchet()` in `phangorn` not find the most parsimonious trees?
I called pratchet() on a dataset with binary values for the characters and it returned a parsimony value of 35. I then calculated the parsimony value for the dataset with a tree that I myself constructed and it was 7. How is this possible at all? I thought that pratchet() was supposed to return the most parsimonious tree on the basis of the data set, but that clearly did not happen in this case. Is there other software that is guaranteed to find the most parsimonious tree?

Please provide a minimal working example of code you have tried. That way a concrete and useful answer can be provided.

Perhaps the tree which you have calculated is not being caught in the search heuristic? Have you tried (as a trial) setting the tree which you calculated as the starting topology, just to check that it is indeed being calculated on by pratchet()?
Alternatively, you could try playing around with the search heuristic settings, like varying from ratchet to NNI, or you could also try using an external program, like PAUP*.

Thank you, I tried these. When I set the start tree to the tree that I created manually, I got the following result:                                                                    Error in edgeMatrix[, k + 1] <- c(p1, e12, e34, p2, 1L) : number of items to replace is not a multiple of replacement length                                                         What does this mean?

I am unsure. Have you verified that you provided a bifurcating newick tree file and that the data matrix has the same number of taxa as the tree?

The number of taxa is the same, but the tree is not fully resolved. So I presume that is then the problem.

For dendrograms with less than about 10 OTUs/taxa/leaves it is possible to evaluate all possible dendrograms using an exhaustive search algorithm. However, for more than about 10 OTUs, heuristics need to be employed to try and "home-in" on and score a smaller set of "better" trees. However, these heuristics are not guaranteed to find the best dendrogram.
Without some specific code it is hard to do more than speculate as to the problem you are encountering, especially where you get an error using data you manually create.


text:
Input normalization in ChIP-seq
If I subtract input counts from ChIP counts (for every gene, since I have one peak per gene) I get negative values for most genes. This makes sense to me, because (as can be seen in the figure) input and IP have same sequencing depth but IP is very biased towards few genes with lots of counts.

Elsevier was informed that the business model which they offered does not meet the demands. How should I implement input normalization?
P.S. The subtraction idea is from Bioconductor

What is your goal?

My goal is to compare RNA-polII ChIP experiments between different stress conditions and different strains

I think if you can clarify your question you will be more likely to get the help you need.

I strongly suggest that you not try to come up with your own package for this when things like CSAW already exist in bioconductor and provide a number of useful normalization options.
For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using bamCompare from deepTools) or alternatively normalize by total coverage (possibly excluding blacklist regions, if applicable for your organism).

package ‘csaw’ is not available (for R version 3.4.3)
I was not very precise. I am evaluating transcription initiation.I am not interested in detecting peaks. I am comparing regions around the TSSs. So far I have been simply normalizing through a spike-in dna from pombe into cerevisiae (my model organism). Any other alternative available before I try to re-invent the wheel?

It's a bioconductor package, so follow the appropriate installation instructions. This is also appropriate for your situation.

If you want to compare conditions (or strains), I don't see why you need Input at all (unless you want to define intervals).
It isn't clear from you question (and subsequent comment) on what you want to achieve, but I can see 2 possible scenarios:

You already have predefined intervals to test on
In your case (Pol2-ChIP), this might be for example gene promoters.
You could compare the Pol2-ChiP libraries directly using standard test for count data as established for RNA-Seq (i.e. edgeR, baySeq etc) in which the intervals are genes (as opposed to promoters or so). No need for input in this scenario, all you need are ChIP libraries for both conditions (don't forget to do replicates).
You don't have predefined intervals
You could defined them through peak calling, in this case you'd in fact need input. A nice approach for this is described in Diaz et al and implemented in deeptools. It poses an elegant method to separate the IP signal from background input, which is probably what you want I suppose. Having defined peaks, you could do scenario 1) using them as intervals.

Case 1 is my scenario, nevertheless.
I have been normalizing the samples with spiked-in foreign dna but I am trying to find out if the inputs make a better job.
I do use limma for downstream analysis, but I try to find the best way to standardize obtaining counts from reads, minimizing the noise.


text:
Is there a database of protein sequences/structure along with their melting temperature?
Such databases have been constructed for example in here or here. But I can't find them anywhere online. Can someone point me to where I can find one available online?
Otherwise I'll have to write to these people and I don't know when/whether I'll get a reply.

ProTherm2 contains thermodynamic properties of proteins and mutants. https://www.iitm.ac.in/bioinfo/ProTherm2/ .
NAR Database issue advertises it as containing melting temperatures.
For the reference - found it with https://biokeanos.com/search?q=protein+%2Bmelting

Thanks! What amuses me is that back when I needed the database, we did contacted the authors of ProTherm for the data and used it. Seems like it's the only available database for this purpose, right now.


text:
Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05
I'm using edgeR for differential analysis. Using glmTreat function I'm detecting differentially expressed genes between Tumor and Normal. I have set the arguments like below:
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=2)
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

This gave me differentially expressed genes with logFC > 2 at FDR < 0.05 and genes with logFC < 2 at FDR < 0.05. [There are few genes with logFC 1.9, 1.89 also]
I have also tried with lfc=log2(2).
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=log2(2))
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

To my surprise I see DE genes with logFC >=1, FDR < 0.05 and also logFC <=1 at FDR < 0.05.
Here is an example of DE genes I also found when I used lfc = log2(2)
GeneSymbols logFC   unshrunk.logFC  logCPM  PValue  FDR
DPYD-IT1    1.14493675  1.560389373 0.67468395  0.002972841 0.018704644
LINC00945   1.144288473 1.525191847 0.693845996 0.00651175  0.038064312
ITPKB-AS1   1.14046752  1.484991704 0.703887505 0.006620218 0.038585862
AC105206.1  1.122495266 1.4871905   0.699167758 0.006545299 0.038194633
LINC02066   1.120660287 1.558458881 0.66016102  0.004344743 0.026424298
AC090985.1  1.10235045  1.58968763  0.64343419  0.002745089 0.017426607
KCNJ6-AS1   1.100399688 1.717028637 0.613324799 0.002224039 0.01438451
LINC01494   1.090270138 1.490899364 0.676937984 0.006770608 0.039348027
AC099541.1  1.080957092 1.53817727  0.651356989 0.008449298 0.048356352
AL133396.2  1.012984207 1.512530231 0.63130643  0.008090259 0.046411928
AC011363.1  1.01207753  1.499708945 0.635740622 0.007226791 0.041797256

I'm very confused with these foldchange cutoffs.
Which one should I use to select differentially expressed genes based on fold change >=2 and FDR < 0.05?

Gordon Smyth has provided a comprehensive answer to your question below, but if any reader is interested in ranking genes by effect size check out TopConfects
(https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1674-7)

Fold-change >= 2 is the same as logFC (log2(fold-change)) >= 1, so your example is doing exactly what you want. logFC is generally easier to think in and work with than fold-change, since (aside from computational reasons) then increases and decreases in expression differ only in sign and are, thus, easier to compare in magnitude (e.g., it's easier to tell that logFCs of 1.5 and -1.5 are the same magnitude change than it is that 3 and 0.333 are the same magnitude change in the fold-change).
For prioritizing the results, typically you would then rank these filtered results by fold-change and relevance (either given background knowledge or for genes with known function).

Thanks a lot for the clear explanation. I guess in my example lfc = 2 is wrong analysis in this case if I need Fold-change >= 2.

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE.
The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset.
There's a few subtleties here:

First, there's the difference between logFC and FC as explained by Devon.

The lfc in glmTreat is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.

The lfc threshold applies to the raw fold-changes (unshrunk.logFC in the output), not to the shrunk versions (logFC in the table).

The intent of glmTreat is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for glmTreat is a bit light, but you can get an idea of how to use it from our workflows.

https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

I understand the simple fold-change cutoff of FC > 2 with FDR < 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline.
For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts.
The glmTreat approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.


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:
Where can I find the DNA sequence of the C. elegans rRNA cluster?
Where can I find the DNA sequence of the C. elegans rRNA cluster?
Ideally, I would be able to do it using the command-line.
I've tried using NCBI but the results are unsatisfactory. I am interested in the DNA sequence that is in the current reference genome of WormBase.
I have tried using the Bioconductor biomaRt package and create a table containing the description but when I grep 'rRNA', I only get proteins doing things with rRNA. I used the following code:
ensembl <- biomaRt::useMart('ensembl', dataset = 'celegans_gene_ensembl')
wormGenes <- biomaRt::getBM(attributes = c('external_gene_name', 'chromosome_name', 'start_position', 'end_position', 'gene_biotype', 'description'), mart = ensembl)
wormGenes[grep('rRNA', wormGenes$description), 6]

[1] "rRNA 2'-O-methyltransferase fibrillarin  [Source:UniProtKB/Swiss-Prot;Acc:Q22053]"
[2] "Pre-rRNA-processing protein pro-1  [Source:UniProtKB/Swiss-Prot;Acc:Q22006]"
[3] "Pre-rRNA-processing protein TSR1 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q19329]"
[4] "rRNA methyltransferase 2, mitochondrial  [Source:UniProtKB/Swiss-Prot;Acc:O62251]"
[5] "Probable rRNA-processing protein EBP2 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q09958]"
[6] "Putative rRNA methyltransferase  [Source:UniProtKB/TrEMBL;Acc:Q5CCJ2]"

The info is actually in the gene_biotype column of my wormGenes data frame. I could select the rRNA genes using dplyr::filter(wormGenes, gene_biotype == 'rRNA'). Unfortunately, no description is available for the genes, so I will have to get it using the gene IDs.


text:
Do any computational phylogenetic methods enable the specification of ancestral states?
Various phylogenetic algorithms estimate ancestral states of a phylogenetic dataset. Is there a way in either maximum parsimony, distance-based methods, or Bayesian inference to indicate what the ancestral states of characters were?

In what respect? Do you mean as in looking at the evolution of a character on a tree? If so, there are algorithms for constrained Brownian motion. If you mean in terms of tree estimation, then knowing the ancestral states implies that you know something about the relationships already - if so, you can use constraint trees for some relationships in ML and Bayesian inference.

Otherwise if you are just looking for ancestral character state estimation, look at the ace() function in the ape package, or make.simmap() in phytools.

Thank you. I had in mind specification of ancestral states in terms of tree estimation, so I will investigate constraint trees in ML and Bayesian inference.

For maximum parsimony, Fitch's method seems simple and intuitive, and can be done in the R phangorn package: Fitch, Walter M. (1971). "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". Systematic Zoology. 20 (4): 406. doi:10.2307/2412116. ISSN 0039-7989. JSTOR 2412116
There are various methodological criticisms of parsimony that you can find in the literature, but its appeal is you can easily understand why the algorithm came to the conclusions it did.
RAxML can do maximum likelihood reconstructions based on a tree and tip values, and I think Mr. Bayes can do a Bayesian version.

Curious, it seems like Systematic Biology has the same paper here

The 70s sure were a crazy time!

Yes this approach is done in parsimony, i.e. you can indicate the outgroup and the ancestral state is calculated from there. I assume this is also performed using 'Beast' because it is a non-reversible matrix. Regular ML and Bayes approaches use reversible mutation matrices therefore the tree is calculated first, you can then specify an outgroup to calculate the ancestral state. However, in a reversible matrix the ancestral state is not used in the calcuation. I hope that makes sense.


text:
Multiway intersection of 8 lists of genes
I have 8 lists of differentially expressed genes for 8 time points (2, 4, 6, 8, 10, 12, 14, 16 hours), the number of genes in each of these lists are not equal.
For example there are 630 common genes between h2 and h4 time points but I have 1143 genes in h2 list and 768 genes in h4 list , so common genes in h2 vs h4 would be 630/1143 and common genes h4 vs h2 would be 630/768.
That would would be great if I have a code to produce a matrix or plot of all possible pairwise combinations of my 8 time points.

cross posted https://www.biostars.org/p/321861/

EDIT: Sorry, I did not see the R tag and wrote the code in bash. I still left the bash example below and added an R version.
Depends on what you want to do, but there are two common ways in which you could store this. Either as wide format:
      h2  h4  h6 ...
Gene1 1   1   0
Gene2 0   1   1
Gene3 1   1   1
...

Where 1 denotes presence and 0 absence, or as long format:
      timepoint
Gene1 h2
Gene1 h4
Gene2 h4
Gene2 h6
Gene3 h2
Gene3 h4
Gene3 h6
...

The best format probably depends on what you want to do. Since you did not specify a method/language, I took the liberty to write a few lines of bash to go from the gene lists to the matrices.
Bash examples assume your files are named h2, h4 etc and are tab separated with gene names in the first column. R examples assume tables is a list of your data frames and timepoints is a vector c("h2", "h4", ...).
Generate table in wide format:

Concatenate gene names of all tables to a separate file
Extract unique values in this file
Perform outer join between the separate file and each table

R:
all_names <- c()
for (tp in timepoints){
  all_names <- append(all_names, tables[[tp]]$genes)
}
all_names <- unique(all_names)

wide_table <- matrix(nrow=length(all_names), ncol=(1+length(timepoints)))
wide_table <- as.data.frame(wide_table)
colnames(wide_table) <- append("genes", timepoints)

for (tp in timepoints){
 tmp <- tables[[tp]]
 tmp$tp <- 1
 tmp <- merge(wide_table, tmp, by="genes", all=T)
 wide_table$tp <- tmp$tp
}

Bash:
Note: you will might need to add headers to this table
for timefile in $yourfiles;do
  # make a file with all gene names
  awk -v tp=$timefile '{print $1}' >> names_file
done

sort names_file | uniq > tmp && mv tmp names_file

for timefile in $yourfiles;do
  join -a1 -a2  -j 1 -o 0,1.2,2.2 -e 0 $names_file <(awk '{print$1,1}' $timefile | sort)
done

Generate table in long format:

For each table, add a column with timepoint
Concatenate all tables

R:
df <- data.frame()
for (tp in timepoints){
  tables[[tp]][[timepoint]] <- tp
  df <- rbind(df, tables[[tp]] )
}

Bash:
   for timefile in $yourfiles;do
      awk -v tp=$timefile '{print $1,timefile}' >> long_table
    done

EDIT2: To follow up on this, here is a simple way (using the package dplyr) to count the number of common genes between all possible combination of time points. Assuming you have your data in long format:
df <- data.frame(genes=c("G1","G1","G2","G3","G3"), tp=c("h2","h4","h2","h4","h6"))

  genes tp
1    G1 h2
2    G1 h4
3    G2 h2
4    G3 h4
5    G3 h6

You can get the co-occurrences using :
library(dplyr)
# Make timepoint ordinal factor
df <- mutate(df, tp = factor(tp, ordered=T))
df <- df %>%
  # Compute all pairwise co-occurrences
  full_join(df, by="genes") %>%
  group_by(tp.x, tp.y) %>%
  summarise(N_genes = length(unique(genes))) %>%
  # Exclude self pairs (gene in h2 and h2)
  filter(tp.x != tp.y) %>%
  # Check for duplicate pairs in reverse order (h2-h4 and h4-h2)
  mutate(comb=ifelse(tp.x < tp.y,
                     yes = paste(tp.x, tp.y),
                     no  = paste(tp.y, tp.x)))
# Remove the duplicate pairs
df <- df[which(!duplicated(df$comb)),]

And visualise it using for example a heatmap:
(ofc it looks stupid because I have a tiny sample where not all combinations are even present)
# Visualize number of genes present in each combination
library(ggplot2)
ggplot(df, aes(x=tp.x, y=tp.y, fill=N_genes)) + geom_tile()

excuse me, i did not quite understand. i have 8 separate lists of genes so how i know about 0 and 1 (absence and presence of a genes among the lists)?

@FereshTeh 0 and 1 were just placeholder values to give the info "geneY is present (1) / absent (0) at timepoint hX" . You could have expression values, ratios, pvalues or whatever variable instead of 1, and NA instead of 0, for example. Does it make sense ?

thanks a lot, actually i just have 8 lists of genes, lets say that each list has one column of genes. i want a pairwise comparison for common genes for any possible combinations of these time points.

Oh sorry, i did not understand the last part about pairwise comparisons. First, if you want to filter common genes, you could just filter rows (genes) which have values available for more than x timepoints. For pairwise comparison, maybe that'd be easier if you provided more details/background. (E.g. do you want to compute a single summary stat per condition and compute it, or is it per gene ?)

thanks a lot, i think per condition. actually i need common genes. for instance i have 630 common genes between h2 and h4 time points, while i have 1143 genes in h2 and 768 genes in h4 itself. so, i have two ways comparison; 630/1143 and 630/768. that would be great if i could do the same for all of my time points. something like having a matrix of percentage or a graph of them.

Oh I think I understand. Could you update your original question with the exact comparisons you want to make (e.g. Compute number of available genes between all possible pairs of combination). Then I can make sure I got it right and update my response.

I updated my solution, please mark it as accepted if it answers your question.

excuse me, i quite did not understand if i have 9 files each with one column of genes, how i can have df (long format). :(:(

I slightly edited my answer. hopefully it is more clear. To generate df, you need to generate the table in long format. I made a section with instructions (bash and R) to generate such table. If you chose R, you might need to edit the code to suite your data. If you chose bash, you just need to read the "long_table" file into R using read.table().

sorry for sure your code will work. i am really sorry but I'm too new in coding. for example for producing long format n r, should i have a df for each time point? i mean 9 df files? i don't know what is tp or should i put my 9 columns in  one file with its own headers for example h2, h4 in header of each column. i really could not figure out. however definitely your code is so practical

No problem with being new to programming. I am afraid I cannot help you anymore than that, but if you get a bit more into programming first, I think it will make sense. There are great free resources out there. Here is a general purpose list:https://github.com/Michael0x2a/curated-programming-resources/blob/master/resources.md
For R I would maybe try: https://www.datacamp.com/
For python in bioinformatics: http://rosalind.info/problems/locations/
For bash: https://www.codecademy.com
Hopefully some of these will be helpful to you :) Best of luck

An UpSet Plot is a useful way of visualising these kind of multi-dimensional overlaps. There's an R package for producing them, with a reasonable set of examples.
The input is a binary matrix indicating set membership:
         h2    h4    h6...
GeneA     0     1     0
GeneB     1     1     0
GeneC     0     1     1
...

And the output is a pair of bar charts, showing the size of various sets, with a dot matrix showing the set membership. An example from the UpSetR documentation:
library(UpSetR)
mutations <- read.csv(system.file("extdata", "mutations.csv",
                      package = "UpSetR"), header=T, sep = ",")
upset(mutations, sets = c("PTEN", "TP53", "EGFR", "PIK3R1", "RB1"),
      sets.bar.color = "#56B4E9", order.by = "freq",
      empty.intersections = "on")

I would simply do this...

# sepcify your timepoints, each timepoint is a list of genes!!
timepoints <- list(h2,h4,h6,h8, h10,h12,h14,h16)
# initialize an empty array that you will fill with your timepoints
vect = vector()
# go over all combinations and fill the array with the set intersection values you want (I guess...)
for (i in timepoints) {
  for (j in timepoints) {
    vect <- c(vect, length(intersect(i,j))/length(i))}}
# reshape your array to a more ordered matrix
your_matrix <- matrix(vect, nrow=length(timepoints))
# assign columns and rows names
colnames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")
rownames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")