text:
Extracting genes from corrplot and adding labels based on high and low corelation
I am doing a gene-gene pairwise correlation, I get a plot and I do see positive and negative correlation
My code:
library(corrplot)
library(Hmisc)
gen = read.csv('NEW_RBP/CORPLOT/TF_NORMAL_NEW_40.txt',row.names = 1)
head(gen)
dim(gen)
gendat <- t(gen)
#gendat <- gen
head(gendat)
dim(gendat)

macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
p_mat <- cor_5$P
row_names <- rownames(M)

pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")

corrplot(M, method = "color",is.corr = FALSE,
         tl.srt = 45,diag = TRUE,tl.pos='n'
         ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dev.off()

How do I extract those gene pairs? One way what I have read and done is:
library(dplyr)
flat_cor_mat <- function(cor_r, cor_p){
    #This function provides a simple formatting of a correlation matrix
    #into a table with 4 columns containing :
    # Column 1 : row names (variable 1 for the correlation test)
    # Column 2 : column names (variable 2 for the correlation test)
    # Column 3 : the correlation coefficients
    # Column 4 : the p-values of the correlations
    library(tidyr)
    library(tibble)
    cor_r <- rownames_to_column(as.data.frame(cor_r), var = "row")
    cor_r <- gather(cor_r, column, cor, -1)
    cor_p <- rownames_to_column(as.data.frame(cor_p), var = "row")
    cor_p <- gather(cor_p, column, p, -1)
    cor_p_matrix <- left_join(cor_r, cor_p, by = c("row", "column"))
    cor_p_matrix
}

#cor_3 <- rcorr(as.matrix(gendat))

my_cor_matrix <- flat_cor_mat(cor_5$r, cor_5$P)
head(my_cor_matrix)

dim(my_cor_matrix)

Delete.na <- function(DF, n=0) {
      DF[rowSums(is.na(DF)) <= n,]
    }

TF_NORMAL_CORR = Delete.na(my_cor_matrix)
dim(RBP_DISEASE_40_CORR)

TF_NORMAL_HIGH = TF_NORMAL_CORR[ TF_NORMAL_CORR$cor > 0.9 & RBP_DISEASE_40_CORR$p < 0.00000001, ]

TF_NORMAL_LOW = TF_NORMAL_CORR[TF_NORMAL_CORR$cor < -0.9 & TF_NORMAL_CORR$p < 0.00000001, ]

Now even if I do that sort of high cut-off filtering I still get considerably high number of correlation pairs like few thousands if positive correlation is considered.
So few of my issues which I would like to get input or help

How do I extract the cluster of regions with positive correlation or negative correlation  ?
If I found certain gene pair/pairs how do I label them ? should I make a vector or map or any other way like if say more than 0.8 correlation I would like to put a label on those regions only
Considering label the gene order is important to be considered, what I have read and found is this how to find the order of the label

tempCor <- cor(gendat)
tempCor <- data.frame(tempCor)
dim(tempCor)
names(tempCor) <- c(1:741)

out <- corrplot(t(tempCor),method = "color",is.corr = FALSE,
                 tl.srt = 45,diag = TRUE,tl.pos='n'
                 ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dimnames(out)

Now with all these how can I extract those patches of high and low correlation pairs and add selected label or highlight them and add label
I m not sure if all these can be done with corrplot! like adding selected label if yes i would be really glad to use that..
My data file data used
"In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples."
Its symmetric so would it matter if i use upper left or bottom right? I guess anyone would do good.
Format - a list for each group if i understand if there is one gene with many positive correlation then that would be considered as group?is that you are saying ?
I have 16 samples in this case but here I using gene-gene. I'm bit curious with a naive question: Is it possible to do gene to sample correlation or if i know what genes are positively or negatively correlated can I trace back without have to plot another heatmap of found genes which are positively or negatively correlated to see their expression in sample ?
I have given the data file it might help you better with my context..

In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples.

I have updated my question with your queries ..hope it helps

Here is a solution using the pheatmap library to cluster and visualise the correlation matrix, then extract the groups from the cluster dendrograms:
gen <- read.csv(file = "TF_NORMAL_NEW_40.txt", header = TRUE, row.names = 1, sep = ",")
gendat <- t(gen)
library(corrplot)
library(Hmisc)
macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
str(M)
p_mat <- cor_5$P
row_names <- rownames(M)

library(pheatmap)
pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")
plot_data_pheatmap <- pheatmap(M, color = rev(macolor), clustering_method = "complete", fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

tree_cut <- cutree(plot_data_pheatmap$tree_row, h = 20)
tc <- data.frame(tip = names(tree_cut), clust_membership = as.character(unname(tree_cut)))
row.names(tc) <- tc$tip
tc <- tc['clust_membership']
head(tc)
str(tc)
class(tc$clust_membership)
pdf("TF_NORMAL_NEW_40_gene_clustcut.pdf",20,15,paper = "a4r")
pheatmap(M, color = rev(macolor), clustering_method = "complete", annotation_row = tc, annotation_col = tc, fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

okay! but how about the question i asked regarding high and low correlation? you saying i should be giving the corelation matrix as input to the pheatmap ?

@krushnachChandra Yes.  Input the correlation matrix to pheatmap.  I have edited my answer.

okay i will give it a try so its the same thing as plotting in pheatmap there won't be any change i suppose..

@krushnachChandra I have updated my answer to use your data and replicate what you have done to explain my method


text:
log rank p value interpretation from gepia 2
How do i interpret the overall survival when the p value reported is not significant going by the standard convention .
Not sure what kind of test it runs while comparing the expression may be one non-parametric test , but the survival is kaplan meier on cox ph model which they have mentioned .
So is it a significant one if i have to tell in context of survival ?

Could you provide a bit of context? From which article did you get this? Why do you think it needs different interpretation when the p-value is above your threshold?

it's from gepia site so im looking for genes in databases which are significant across patient sample ,boxplot shows significant but the survival is not significant if i understand ,why is it so?

The plot on the left shows that the gene has very different expression between two groups. However, according to the graph on the right, these high expression is not directly related to survival because both groups have similar survival time.
The gene expression so high might be to compensate for other deficiencies or problems.

yes both of them have similar survival  which is reflected in the log rank value ? in other words if i see the log rank value is not significant ,then i shouldn't consider it..

@krushnach You shouldn't consider it for what? Please [edit] the question to include more details (I don't even know which gene is: LAML or zbtb16??), perhaps a link to the result on the website is possible

okay i will do it. i cannot link it as it the result is stored and rendered there but here is the link http://gepia2.cancer-pku.cn/#survival the data set im using LAML gene is ZBTB16 you can have a look..

Ok, that's for the survival curve and the other plot were does it come from?

http://gepia2.cancer-pku.cn/#analysis this one

Ok, know I know how did you get these plots. But why do you need them ? Why are you looking at these plots? What are you looking for

as they were getting differentially expressed in leukemic stem cells and blast cells so was looking in TGCA how the survival is affected ..

source: Stack Exchange

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

sort: votes


text:
List of human transcription factors
I'm looking for a list of human transcription factor (TF) genes.
Where can I find a TF database or how can I retrieve this list from other databases?

Lambert et al. (2018, Cell. 172(4):650-665) provide a manually curated list of human TFs at this website.

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


text:
Receptor-ligand database
How can I find ligands of receptors (or vice versa)?
For example, I want to find the receptors of IL-8 (CXCL8).

CellPhoneDB is a publicly available repository of curated receptors, ligands and their interactions.
The interactions can be searched or downloaded.

Have you tried STRING?
It's a publicly available protein-protein interaction database, so it won't be a complete picture.
Ingenuity Pathways Analysis (IPA) is really nice but is a paid service. Ask your institution if they have a license for one.

Can you include a link to IPA please.

https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/

CellTalkDB is another option for a searchable ligand-receptor database. According to the webpage it is curated from other databases and the literature. Details are available in the publication.


text:
Normalize before using adonis?
I'm using adonis, an implementation of MANOVA in R (vegan package) to understand how continuous (temperature, pH) and categorical (Site ID, Month) variables affect my 16S dataset. Applying it using my metadata 'as-is', I see that ~80% of the variation (R2) is significantly explained by one variable - Site ID.
However, after suggestions from colleagues to normalize my data (as shown here for 'min-max' normalization), the same test not only shows Site ID as the third most important variable (in explaining the most variation significantly), it only shows ~31% variation being explained from the first three variables.
Any opinions on best practices for this kind of analysis?

Hi and welcome to the site. I'm afraid cross posting is very much frowned upon on the SE network. Please either delete your question from one of the two sites, or edit it so that it is specifically addressed to the subject matter of each site. For instance, for us, we would need to know what kind of data you're dealing with ("16S dataset" isn't very specific).

Also what kind of normalization did you do?  What test are you using? t-test or some package normalization and test procedure?

@Llopis that would be a "min-max" normalization, as what's found here: https://vitalflux.com/data-science-scale-normalize-numeric-data-using-r/
adonis is documented here: http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/adonis.html
It uses something directly relatable to a MANOVA.

@AndréSoares Could you edit that comment into your question? Also if you could add what kind of data are you using (like terdon asked)? Why did you choose the "min-max" normalization?

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


text:
Assembling sequence data generated by RADseq
I'm working on a project where I have to assemble sequences generated by RADseq. At the end I hope to compare two species of woodpeckers in Sri Lanka by using SNPs.
I tried to assemble it using Stacks pipeline, but it crushed on the way with:

*** Error in `ustacks': malloc(): memory corruption: 0x000055db58024be0 ***

Are there are any other tools or pipelines that I can use for assembling my data instead?

There are plenty, see a list here.
If you want a better answer, write a specific question, write what you tried, why Stacks doesn't work and what is the desired output.

@ShanikaLasanthi I tried to make your question clearer, it would be good to add the details about how you run stacks and what kind of error message you got before it crushed.

I don't think you should be assembling RADseq, you should be aligning it.

I am getting this error saying "*** Error in `ustacks': malloc(): memory corruption: 0x000055db58024be0 ***"
what am I supposed to do now?

@ShanikaLasanthi Which command did you use?

I used the code denovo_map.pl -s trim_ln7_brds_KD03S01.fastq -o ./stacks -b 1 -B trial1_radtags -t -T 15 -m 3 -M 2 -n 1.

I figured out what happened. My sequences were not in the same length,therefore it interfered with the algorithm giving this error. once I trimmed the sequences into the same length, pipeline was implemented successfully without giving any error.
Thank you all for the suggestions and the advice.

@ShanikaLasanthi it would be nice if you write here an actual answer for your question. It seems that the solution is not very intuitive and it could be very helpful for others to see the answer.

Well, first my sequences were in different lengths between 50-100bp. I got the error when I used these sequences in the pipeline. When I checked the denovo log file for the run ( you can find it in the output folder you have mentioned for the run) there was a warning saying "different sequence lengths detected, this will interfere with Stacks algorithms."
Therefore I trimmed the samples into 80bp and filtered out sequences shorter than 80bp.
After that I did not get the error message and the pipeline was implemented successfully.
And I hope this helps.

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


text:
How to obtain clusters of hierarchical heatmap when using Python?
Is there a good way of obtaining the labels (e.g. genes) within individual clusters that haven been clustered hierarchically in Python (preferentially, but not necessarily, by seaborn)?
I found these difficulties highlighted elsewhere (here, and here), and would hope for some simple option (e.g.: like in MATLAB where cutoff can be specified, and corresponding genes would be returned), without potential risk that the plotted hierarchical heatmap differed from the cluster assignment ( see here)

This question needs more information; it's not specific enough. Python is a programming language, not a toolkit. At the moment, it's similar to asking "Is there a good way to find the names of common taxi companies that are found in cities in Europe (preferentially, but not necessarily, French cities)?"

@gringer; either solution or library would be fine as long as it can be called from Python and does the above job and works well (which would exclude for me solutions that go through Python executing compiled local code through bash, or calling local instance of R or MATLAB); perhaps: "is there at least one country in Europe, where it is possible to order a taxi (preferentially, but not necessarily, by phoning the taxi company)?"

Clustering like this is typically done with scipy. Here's the code we use in deepTools (original context here):
from scipy.cluster.hierarchy import fcluster, linkage
Z = linkage(some_matrix)  # You might want to set `method` and `metric`
groups = fcluster(Z, nGroups)  # You might want to set `criterion`

groups is then a vector containing the group assignment of the input matrix according to the clustering. As an aside, this is exactly the same as what's shown in the last link in your post, so I'm not sure how you didn't find it.
BTW, at least in python, hierarchical clustering is painfully slow with large datasets.

thanks; I had indeed been wondering about Python's usual libraries for clustering, including scipy, being awfully slow for hierarchical clustering


text:
Metagenomic shotgun data with internal control
For 24 human stool samples I have metagenomic shotgun reads from an Illumina platform. The reads were filtered and mapped against a bacterial species library and specific maps to species were kept and counted. Importantly, before DNA extraction from any sample it was spiked with a solution of a fixed DNA concentration of Thermus thermophilus, a species not occurring in the human microbiota. The reads from T. thermophilus allow me to convert relative abundance estimates back to abundance in the stool sample, to absolute abundance.
I want to test whether T. thermophilus is indeed a good internal control that removes confounding by converting to absolute abundance. To test this, I compare between-species correlations on the relative abundance scale with their correlations on absolute abundance scale. I would expect less correlation in the latter. I first take the raw read counts of a random set of species, I divide counts by a sample-specific size factor estimate (total number of reads r in sample i / geometric mean of r's) to adjust for library size, then I do arcsinh transform (to stabilise variance) and then I look at between species correlation. Next I take again the raw read counts, correct by the size factor, then correct for an estimate of the coverage ( estimated proportion of T. thermophilus in sample i ), and I do the acsinh transform. Again I look at between-species correlations.
I compare species-pair correlation before and after T. thermophilus correction in a scatterplot (see attached) to visualise the effect of the correction. I expected that the point cloud indicates a slope <1 reflecting that previously strong correlations have weakened on the absolute scale. But that is not what my result shows. Indeed, some weakly positive correlations have become stronger. How does this make sense? What other ways are there to test validity of an internal control?

The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control.
As you can see on figure 3 of the linked article the correlation on the absolute scale tend to be higher. So I would say this is within expectations.
As other way for internal controls are blanks (Just water to see if there are any contaminants), and technical replications, to make sure that the measure is correct not some problem with picking some amplifications instead of others...
Also you don't mention it but the tool and library used to impute the microorganism is also important. You can check with other tools to see if the estimation is similar.


text:
DMPfold: number of models produced using "10 500" feature in "run_dmpfold.sh"?
I ran this command:
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./folder 10 500

From my previous discussion, I understand that I should run with default "run_dmpfold.sh" parameters i.e. 3 50; however, since I have already run the 10 500 I have some questions.
I see that after generating "sequence_1.pdb to sequence_500.pdb" files, these files are no more in the folder and, now I see "ensemble_1.pdb".
Will there be more ensemble files like "ensemble_500.pdb" and "sequence_n.pdb" files? I mean when can I expect "final_1.pdb' file? Also, I see files like "CLUSTERS_001.pdb".

The sequence_n.pdb files are temporary and are combined into an ensemble.i.pdb file after each iteration i. That file has many structures in. After the run is complete there will be a final_1.pdb file. Don't worry about the cluster files. If you are interested in what they do, consult the run_dmpfold.sh script and see what it is doing line by line.
As I have said multiple times elsewhere, you should run the script with default parameters and let it run to completion. Then use the final_1.pdb and any other final_*.pdb files as the output. If you want more structures, use the ensemble.i.pdb files.


text:
1000 genomes missing SNPs on chr X in GRCh38
I am trying to get a reference together for imputation using the 1000 genomes data set. The GRCh38 version of the data set was made available earlier this year, but seems to be missing a significant number of SNPs on the X chromosome (~96%).
Is there any resource that can explain why? I have tried searching on their home page with no luck.
The data sets in question are GRCh38, GRCh37.

Please do not cross post to BioStars, Stack Exchange and 1000 Genomes helpdesk as it can create duplicate work:
https://www.biostars.org/p/399254/

As I told you on 1000 genomes helpdesk:
Because of the ploidy settings we used for our variant calling, we only called variants in the the pseudoautosomal region on the X chromosome. We are working on settings to correctly analyse X variants in both males and females and plan to analyse X properly soon.

I think these are biallelic snps:

The set is restricted to biallelic SNVs and INDELs.

If you check the liftover version, the chromosome X is complete...


text:
Biopython Genbank writer not splitting long lines
I am parsing a csv file of annotated sequences and using Biopython to generate Genbank files for each. I want to add annotations of the sequence features.
My output file shows features listed without the correct line breaks. Other software is then unable to parse the names of the features
misc_feature    58..60  /note="Oxidation"
misc_feature    100..102  /note="Oxidation"
misc_feature    241..243  /note="Oxidation"

It should be formatted like this
misc_feature    58..60
                /note="Oxidation"
misc_feature    100..102
                /note="Oxidation"
misc_feature    241..243
                /note="Oxidation"

My section of code adds a space between the location and the /note="" but the Biopython writer doesn't add line breaks.
for i in range(len(start_list)):
    feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
    feature.key = "misc_feature"
    liab_type = "{}".format(liability).split(":")[1]
#need to add 1 to start list since csv counting is off
    feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i]) + "  " + '/note=\"{}\"'.format(liab_type)
    container.features.append(feature)
    i +=1

The Biopython source is toggling the wrapping of the feature qualifiers on no_space_key and Bio.GenBank._BaseGenBankConsumer.remove_space_keys
How do I enable wrapping of text to format output correctly?
    def __str__(self):
    """Return feature as a GenBank format string."""
    output = Record.INTERNAL_FEATURE_FORMAT % self.key
    output += _wrapped_genbank(self.location, Record.GB_FEATURE_INDENT,
                               split_char=',')
    for qualifier in self.qualifiers:
        output += " " * Record.GB_FEATURE_INDENT

        # determine whether we can wrap on spaces
        space_wrap = 1
        for no_space_key in \
                Bio.GenBank._BaseGenBankConsumer.remove_space_keys:
            if no_space_key in qualifier.key:
                space_wrap = 0

        output += _wrapped_genbank(qualifier.key + qualifier.value,
                                   Record.GB_FEATURE_INDENT, space_wrap)
    return output

Glad you found the solution!, The way to signal that the question is solved is posting your own question and selecting it as the one that resulted more helpful. This will unify the style for all questions and answers and make it easier for people to find the solution to your question.

I read the Biopython Genbank record source code carefully and realized I wasn’t using the Qualifier class correctly.
After reading the source code again, I realized I was using the data structure incorrectly.  I need to use the Qualifier class to add a key and feature.  This is formatted correctly and solves the problem.
feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
feature.key = "misc_feature"
#need to add 1 to start list since csv counting is off
feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i])
liability_qualifier = Qualifier()
liability_qualifier.key = '/note="'
liability_qualifier.value = '\"{}\"'.format(liability).split(":")[1]
feature.qualifiers.append(liability_qualifier)

Thanks for posting back and explaining the answer to the rest of us


text:
How to convert a list of genes to one element per row?
I got two lists of about 300 genes from my GO analysis, like the following:
Mtss1/Sox11/Fbn2/Mmp2/Tcf7l2/Six2/Ahr/Fat4/Lef1/Stc1/Acvr2a/Glis1/Nog/Igsf10/Twist2/Igf1/...

How can I convert the list to one element per row (for Venn diagram  on Venny 2.1), as:
Mtss1
Sox11
Fbn2
Mmp2
...

$ echo "Mtss1/Sox11/Fbn2/Mmp2/..." | tr '/' '\n' > list.txt

For more information on tr: https://en.wikipedia.org/wiki/Tr_(Unix)

It works well in my terminal, thank you so much!

In R you can use tidyr's function separate_rows.
install.packages("tidyr") ## if not already installed
library(tidyr)
GO_results<-read.table("genes.txt")
separate_rows(GO_results, GO_results$V1, sep = "/")

Many ways, assuming you're on a *nix system.

tr. Alex already gave you a fine tr solution, but here's another approach assuming you have an input file:
tr / '\n' < golist.txt > newlist.txt

sed
sed 's|/|\n|g' golist.txt > newlist.txt

Or, to edit the original file:
sed -i 's|/|\n|g' golist.txt

Perl
perl -pe 's|/|\n|g' golist.txt > newlist.txt

Or, to edit the original file:
perl -i -pe 's|/|\n|g' golist.txt

grep (GNU grep, the default on Linux systems)
grep -Eo '[^/]+' golist.txt > newlist.txt

Thank you all very much!

@Lin you're very welcome. If one of the answers solved your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. Since Alex answered first, you might want to accept his.


text:
Understanding Single Cell RNAseq Plots
I have just started some single cell RNAseq analysis with 500 cells. I have used Seurat for analysis. Following the guideline I have generated few plots. Now I am trying to understand how to read those plots. Can the community provide their input in the following plots. To give a brief background of the experiment. This data was processed in MinION. Its from human brain sample from our collaborator. In their shortseq data they identified 15 clusters. However, both the libraries were processed independently. The main question I am trying to answer is celltype specific isoforms. but in my dataset I can find only 6 Cell clusters. So, I am trying to understand why I am getting such a small cluster? Is it because I stared with only 500 Cells? Also I want to know what these two PCA plot is indicating. Thanks

What are your specific questions?

@haci am trying to understand whether the data looks good or not. This data is from long read RNAseq and I had never worked on even short read scRNA seq.

The data looking good or not is not a specific question, it is subjective. For example how many cell types or states do you expect and does it make sense (regarding your research question) to end up with 7 clusters in the UMAP plot? Unless you provide some info on your experimental set-up and research questions you cannot get useful feedback.

@haci Thanks for your input. I am editing my question

How many PCs did you use? And how many cells are there in the dataset where 15 clusters were identified?

Identifying 15 clusters from 500 cells is challenging. It will depend on the relative abundance of each cell type within the tissue. You can use the cell type frequencies found in the shortSeq data set to see what to expect from a 500 cell dataset using the howmanycells tool
About the pc plots:
The first is an elbowplot which explains how many variation is explained by each PC. The knee basically defines when components become less informative.
EDIT: I noticed that the y-axis of the elbowplot starts at 4 suggesting that the dimensionality of your data is most likely ~10 dimensions. But I expect many of the clusters that don't show up on your Umap plot to be present at very low frequencies.
The second plot is a jackstraw plot and provides a visualization of a statistical test that calculates the significance of your principal components. This one shows that all your components are significant. You can combine this with the elbowplot to estimate the dimensionality (10-15 dimensions) of your data and decide which components to use for downstream analysis. In general the more they dropoff to the dotted line the less significant they become.
While these plots are useful as a first impression, always rely on biological information and common sense to decide what information (here: components) to include in downstream analysis such as Umap, clustering, etc. For example: usually you don't want to include components mainly defined by ribosomal, mitochondrial or hsp genes.

Whats the importance of second plot? What can I read from this plot?

@user3377241 updated my comment

The first two plots are described on the Seurat PBMC tutorial. Specifically:

‘Elbow plot’: a ranking of principle components based on the
  percentage of variance explained by each one (ElbowPlot function).

Based on the tutorial, you should use at least 10 PCs.

The JackStrawPlot function provides a visualization tool for comparing
  the distribution of p-values for each PC with a uniform distribution
  (dashed line). ‘Significant’ PCs will show a strong enrichment of
  features with low p-values (solid curve above the dashed line).

Based on the tutorial, you should use more than 15 PCs.
No one can tell you the right number of clusters. If you already identified 15 clusters in a different dataset, then you should have cluster markers for those genes. You can check those genes in your dataset.

The first two plots are used in order to estimate the number of PCs to be used in later stages, for example clustering (this was split in to two functions in Seurat v3). The number of PCs selected would have an impact on the number of clusters obtained. The more PCs the more information for downstream applications. It would be a good idea to go over the genes within each PC, a trained eye can tell what kind of biology (or technical artefacts) contribute to the PCs.
Moreover, FindClusters() has a resolution parameter, which enables the number of communities (clusters) that is generated. In Seurat v3, there is a wide selection of different clustering algorithms but recommending one over the other would be beyond my expertise.
If you have a comparable dataset generated with some other technology/method, you can give a try with "integrating" the two datasets. An alternative non-Seurat approach that aims to do the same is MNNCorrect from the scran package.


text:
How to retrieve the subcellular location info from uniprot?
I have a list of proteins with their gene names/symbols. I would like to identify which of them are extracellular proteins or plasma membrane proteins.
I can get this information in UniprotKB by searching the proteins one by one. Because I have hundreds of proteins, I wonder whether I can do a systematic search of the UniprotKB in RStudio. I tried the UniProt.ws select method. But I couldn't even get the example code in the package work. I downloaded the uniprot_sprot.dat file. But I don't know how to read this file correctly in RStudio. I used tab-delimited format. It turned out there is only one variable in the file.
Can anybody tell me what I should do to solve the problem? I am pretty new to any computer languages.

Can you give us an example ID and where on the UniProt website you see the field you're looking for? Examples always help make the question more specific and grounded.

Welcome to the site. I edited a bit the question to make it easier to find it and read. You could use gene ontology cellular component to know where a protein can end. Also, show how did you try to read the file and why did you download that file? I couldn't find anywhere that it has the information you seek.

Since you are already on UniProtKB, the simplest way will be if you add this information to the output table and download it as a that you can import and filter in R.
To get additional annotations, click the button labelled "Columns" on top of the UniPortKB website.

This allows you to add a wide variety of annotations such as subcellular location to the output. Tick the checkmark for subcellular location and save on the top right.

Then you can download the proteome-wide annotation for example as a tab-separated list using the "Download" dialog next to the "Columns" selection. This list can be imported into R.

Hi PPK,Thanks for the information on how to download the subcellular location info from UniProtKB. The file seems very large to be downloaded. I am in the process of downloading.

After my post of the question, I went to try the UniProt.ws package in the RStudio again. I realized that there is a mistake in the example code in the package. The up object is from a human database, but the keys the example provides are the ENTRZID from mouse proteins. When I tried human ENTRZID numbers, it worked. The only problem for the select method in the UniProt.ws is that it can't retrieve the protein information for a large batch of proteins. The maximum number that worked for me is about 80 proteins per time.


text:
Mapping GPL Platform Number to Microarray Name
I am trying to download a dataset from GEO using the GEOquery R package. Using this package, I can retrieve the GPL platform number(s) for the dataset, and I hope to use a CDF file from BrainArray on the dataset in order to undo any existing normalization performed, so that I can have multiple datasets that are normalized in a standard manner.
In order to do this, I need to know what the microarray platform name is. This information is usually provided on the website online (for example, here, you can see that next to GPL96, it says HG-U133A, which is the CDF I would use from BrainArray).
However, I am wondering if there is any way to parse this platform name using the GEOquery package in R itself rather than having to manually check for every dataset and compare it to the BrainArray website?
If it is not possible to do this in the GEOquery package, is there any online tool for mapping GPL to platform name? I found this online which has the mapping for a decent amount of GPL numbers, but not all of them, meaning it sometimes crashes. I was wondering if there is a foolproof way to do what I am trying to do.
Thanks for the help in advance!
Edit:
I know I can access the highlighted portion in the first image below (the GPL platform number) through GEOquery. I'm wondering if I can access the name of the platform (highlighted in the second image). If not, my alternative option will probably be to scrape the web page to find this information, but if its possible to do so through GEOquery I'd rather do that.


text:
How to remove degenerate sequence (N4) at the 5' and 3' of the reads?
I have to analyse GSE114327. For that, I use GSE instruction based on GSM3139597.
In that data processing section has written:

FASTQ reads were subjected to adapter removal (cutadapt) per manufacturer's instructions (NEXTflex Small RNA Sequencing Kit v3).

So, I use cutadapt for trimming. Then

The degenerate sequences (N4) at the 5' and 3' end of the reads were removed

based on the instruction. But I don’t understand that sentence and I don’t know how I can remove degenerate sequence (N4) at the 5' and 3' of the reads.

From the link you posted it looks as though they used UMI tools as per the section which says
'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'.  The link above has instructions for installation and use.
From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read.
Looking up the protocol for the nextflex kit used (See Link) confirms this:

Data Analysis
  The 3' and 5' adapters included in this kit both contain 4 random bases that will appear immediately 5' and 3' to the insert in sequencing data. The presence of these random bases should be considered when choosing an alignment strategy. When using “end-to-end” alignment, we
  recommend processing data in the following manner:
  1. Clip the 3' adapter sequence (TGGAATTCTCGGGTGCCAAGG).
  2. Trim the first and last 4 bases from the adapter-clipped reads.
  3. Perform alignments as normal.

It seems that N4 is a 4 base random barcode at both the 5'& 3' end of the read - these are cut from the read and this metadata is added to the read name so that every read can be uniquely identified.
Edit:
You can install UMI tools via:
conda install -c bioconda -c conda-forge umi_tools
or
pip install umi_tools
Usage instructions can be found at:
https://umi-tools.readthedocs.io/en/latest/reference/extract.html

Usage:
  For single ended reads, the following reads from stdin and outputs to stdout:
umi_tools extract --extract-method=string --bc-pattern=[PATTERN] -L extract.log [OPTIONS]
For paired end reads, the following reads end one from stdin and end two from FASTQIN and outputs end one to stdin and end two to FASTQOUT:
umi_tools extract --extract-method=string --bc-pattern=[PATTERN] --bc-pattern2=[PATTERN] --read2-in=[FASTQIN] --read2-out=[FASTQOUT] -L extract.log [OPTIONS]

Dear@Graeme, let me ask that how can I find [PATTERN] for GSM3139597?


text:
rentrez to retrieve NIH files
I'm trying to use rentrez package to retrieve NIH files, and ran the following code (I followed what was on the rentrez tutorial):
pubmed_search <- entrez_search(db = "pubmed",
                    term = "Case Reports[Filter] AND cardiovascular disease AND English[lang] AND 2009:2019[PDat])",
                    retmax = 792711,
                    use_history = TRUE)

for(i in seq(1,10000,100)){
    recs <- entrez_summary(db="pubmed",
                           web_history=pubmed_search$web_history,
                           retmax=100,
                           retstart=i)   cat(
seq_start+99, "sequences downloaded\r") }

However, instead of the 10,000 texts I expected, I only downloaded 100 files. Could so

Welcome to Bioinformatics.SE! It looks like you posted the question before you were done typing. No worries, just use the edit button to update it. Also, please indent code with 4 spaces or wrap it in triple backticks (\```) so that it displays in a legible way.

Also, include an API code in your Pubmed search script to increase the query limit.

I don't program in R, but replace the last line of you code with the following and I suspect you will find it works,
for(i in seq(1,10000,1)){ recs <- entrez_summary(db="pubmed", web_history=pubmed_search$web_history, retmax=100, retstart=i) cat(seq_start+99, "sequences downloaded\r") }

Rationale I would guess seq is similar to the linspace function in Python, and put simply 10000 (total sample size) / 100 (the step size) = 100, 100 being the number of downloads you observed. In the code above its 10000/1. I would also suspect "step size = 1" is default, if this was correct then seq(1,10000) would give you your required output.
As a random aside I do like the presence of curly brackets defining the loop, it reminds me of Perl.


text:
Compare multiple alignment results' aligned bases
I have aligned a nanopore data set to a reference genome with graphmap, minimap2 and BLASR. The alignment results are stored in BAM files.
I would like to do some concordance assessment, looking at the number of base pairs that are mapped to the same position by different aligner.
Are there any tools or ways that allow me to do this effectively?

is this a DNA or RNA dataset you aligned to the reference genome?

it is the RNA dataset from https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md

I'm not aware of any tools that do this out-of-the-box. Assuming some programming experience, it would be fairly straightforward to implement this in Python using pysam (or htslib bindings for your favorite scripting language). In particular, the pileup is a very useful data structure.
Ultimately, you'll need to store an index of reference sequence positions by read name and position.
read 1, offset 0 --> chrom 3, offset 248879
read 1, offset 1 --> chrom 3, offset 248880
read 1, offset 2 --> chrom 3, offset 248882
...
...
...

Once you've done this for a pair of BAM files, computing the concordance statistics involves trivial set operations.
There are some important practical considerations though. If you implement the approach completely in memory, it will only be tractable for small-ish BAM files. It wouldn't take very many reads (relatively speaking) to overwhelm available RAM on a laptop or desktop. For large BAM files, you might even overwhelm the available RAM on a big-memory server. Alternatively, you could store the indexes in an on-disk database such as sqlite3. This would keep RAM requirements low, but would slow things down substantially.
I think the most reasonable approach would be to process the BAM files in batches: grab all the reads in a 100kb region (or whatever size makes sense), compute the index, store the congruence statistics, move on to the next 100kb region and repeat, and then aggregate statistics when all regions have been processed.

I don't know of a single tool that can read multiple sets of BAM files and compare them. But you can use the Sam_Reader python3 module to create a csv file that summarizes each position in the genome for each of your 3 results and then compare them manually.
It should be installable with:
pip3 install sam_reader

And then generate the CSV with
from sam_reader import Sam_Reader
my_files = Sam_Reader('my_bam_file_or_folder/')
my_files.per_base_stats(write_file='new_file.csv')

*but FYI this module requires pysam, which does not work on windows

You can use the bam_alignments_compare.py script from the wub Python package.

Welcome to Bioinformatics.SE, and thanks for your answer. Note that when linking to external sites, please provide some context. Internet links tend to have a short half life, so describing, quoting, or embedding the most relevant content from the link helps preserve the information. See the Provide Context for Links section of this help page for more details.


text:
ChrQ definition
What is ChrQ?
I assume it then means chromosme Q or something like that yet when I google for it I get very weird answers, but as I saw that a couple times when looking on freebayes pipelines manuals I wanted to be sure.

Please [edit] your question and include an example (or a few) of where you've seen this. It's almost certainly just a place holder, as Devon suggests, but we could give a more definitive answer if you could link us to the places you've seen this used.

Answer from @devon-ryan and @terdon, converted from comments:
Are you sure it's not just a made up chromosome name that's useful in examples since people won't blindly copy/paste the command and have it work?
It's almost certainly just a place holder, as Devon suggests, but we could give a more definitive answer if you could link us to the places you've seen this used. If you only see it in relation to freebayes then yes the author is using it as “foo”.
[please update / edit to improve this answer if you can]


text:
How to interprete PlasmidSeeker results concerning plasmids?
I've been interested in plasmid detection in bacterial genomes. I tried plasmidfinder and this software found only one plasmid in each of my two bacterial genomes. I was dubious about these results, because I read many genome announcements which found something like from 6 to 10 plasmids per genome around 5 -6 MB. I then used PlasmidSeeker and what I received in the result file was a list of 9 or 11 so called "plasmid clusters", each containing a few plasmids with complete ID. I know that my genome cannot contain tens of plasmids, so my question is which one(s) of them is (are) the real one(s)?
By the way, can you suggest other bioinformatics tools specially designed to detect plasmids in bacterial genomes?


text:
How to predict protein binding to ligands?
I want to predict from protein sequence if the protein binds to metal, nuclear or small ligand. How can I do this ?  Which features are relevant if I want to use them in a machine learning algorithm ?

Welcome to the site! Could you please explain what features/data do you have available? Also it could be helpful too if you provided which machine learning algorithm do you want to try. Also could you clarify what does it mean "metal ligand" ? I've found this from wikipedia, but I doubt that this is what you mean

If you are training a machine learning algorithm you probably want to train against things you know to be true, rather than making predictions from the sequence and training against the predictions.
In this case you will want to use protein structures with the relevant ligands bound as your data. If you go to the RCSB PDB advanced search you can search for entries with any ligand present, specific ligands present, protein and DNA present, etc. Then you can extract the FASTA sequences for these entries.
See also databases like MetalPDB for specific ligand types, and one of our papers for an example of creating a dataset of metal binding proteins.
Features to use would be the raw sequence as a one-hot encoding on amino acid types. You could add other sequence level features such as secondary structure predictions, disorder predictions or physico-chemical properties. However a decent network would be able to predict much of this if it needs to. You could create alignments and put in a sequence profile or covariation data, which would be giving the network more information.


text:
The visualisation of a list of genes on URD object
Seurat R package has some functions like FeaturePlot, DimPlot and DoHeatmap by which we can plot the expression of a list of genes on cell clusters. I am working with URD that likely does not have such options or I can not find that. URD has plotDot but would be clear only for a few number of genes (please look at the picture). Could somebody please help me in plotting the expression pattern of a list of genes in URD? As I am not good in R at all to write my own function or adopt Seurat functions to URD

For example Seurat DoHeatmap gives such a picture by which we can see these set of genes are being expressed on which group of cells
DoHeatmap(object ,genes.use=a vector of genes)

Or we can plot a set of genes on seurat feature lot to they are being expressed more on which cluster
FeaturePlot(object = pbmc_small, features.plot = 'PC1')

In URD only by
plotDot(object.6s.mnn, genes = a vector of genes, clustering="Infomap-60")

We can plot a set of genes but as you are being in the first picture that would not be clear much
My question is how I can plot a vector of genes on URD object
for example Seurat DimPlot function can visualise the expression pattern of a group of genes on cell clusters but when I am trying to do the same with URD, returns this error
> plotDim(axial.6somite, c("DDB_G0291149","DDB_G0291029"), plot.title="SHHA (Notochord marker)")
Error in names(data) <- cells.use :
  invalid to use names()<- on an S4 object of class 'dgCMatrix'
In addition: Warning messages:
1: In if (label.type == "meta" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
2: In if (label.type == "group" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
3: In if (label.type == "sig" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
4: In if (label.type == "gene" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
>
>

Actually only one gene can be visualised by this function on URD

Which plot do you want? Could you post an example from the similar function from Seurat? You say you are not good at R, did you try learning with a bootcamp or with online books (such as this one?

Please [edit] your question and give us a minimal, complete, verifiable example that we can use on our machines to reproduce your code. Those users who can answer will need to have access to a test dataset they can use to test their answers. Not your entire data, just a small part of it, or a similar dataset that we can play with.

Sorry, I know something in R but not as good as to be able write a function in a deadline

I just added URD clustering labels for each cel to seurat object and then I did DoHeatmap by using group.by="URD clustering labels"


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:
How to set new identities for new clusters based on a gene expression in Seurat?
I have a group of cells that have yellow (EYFP) expression and cells that don't. I want to define the yellow and non-yellow genes in two cell identities.
I clustered the cells and then created Violin Plot for EYFP expression. Then I did FindMarkers between clusters that express high EYFP with low EYFP expression. However, this is not a very accurate way to isolate cells of high and low EYFP expression. I need to use ident to separate EYFP +ve and EYFP -ve cells (but don't know the exact syntax). Then FindMarkers between these two new identities.
I need to subset a certain set of barcodes that express a marker gene (in this case Olfm4).

Welcome to Bioinformatics! Could you please clarify what have you tried? This will help us to give a better answer to your question. Also, it seems like the text of the question is incomplete, what do you want to compare the yellow and non-yellow genes or something else?

you posted another question that looks like this one but completed. You can [edit] your questions to improve them (and you can delete this question if you do not longer want it).

I subseted a group of cells and then created a Seurat object out of the subsetted cells. I clustered the cells and then created Violin Plot for EYFP expression. Then I did FindMarkers between clusters that express high EYFP with low EYFP expression. However this is not a very accurate way to isolate cells of high and low EYFP expression. I need to use ident to separate EYFP +ve and EYFP -ve cells (but don't know the exact syntax). Then FindMarkers between these two new identities.

Sorry, but can you [edit] your question to include this? What is your question? I don't understand what you are trying to do or what are you having problems with. If you need help with the syntax, please post your code (to know which specific problem you have with the syntax).

I need to subset a certain set of barcodes that express a marker gene (in this case Olfm4).

I need to subset a certain set of barcodes that express a marker gene. Olfm4.fetched<-FetchData(KRAS,vars='Olfm4') to fetch the barcodes and then to create a Seurat object I use Olfm4<-CraeteSeuratObject(Olfm4.fetched). But then I get errors:  Olfm4<-CreateSeuratObject(Olfm4.fetched)

You have to define your ident it has to have the structure of a table. You put the cell names as names and clusters as value (for instance cell 1 = YEPF+ cell 2 = YEFP-):
pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

Then you can run the FindAllMarkers function on this object:
FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)


text:
Protein SCA and DCA
Could you please suggest suitable R packages for statistical coupling and direct coupling analysis of protein sequences?
I know there are Python packages for those: SCA and DCA.
I was looking for similar packages in R so that I can continue with my present R workflow.

I think R's official blog would be helpful for this

Welcome to the site Tamoghna. Could you please explain what are your reference packages in Python? Also, which resources have you looked in? Did you find anything helpful in Bioconductor or similar sites?

SCA: http://reynoldsk.github.io/pySCA/
DCA: https://github.com/giopina/pydca

@TamoghnaDas please [edit] the question to add these references. (someone can overlook them).

Its an interesting question, because as you say this normally the area of Python. The context of the question would be useful, convolution could be an issue

While my suggestion doesn't technically fit your qualifications, I do think it may be useful. I bring this up because you mention you want to continue with your R workflow. I am assuming you don't have anything against Python, but just want to stick with R to keep a tidy project.
The R package reticulate (homepage) allows for a very easy interface between Python and R. I use it to run the Python package scrublet to remove doublets from scRNA-Seq data before analysis with a popular R package.
Just for reference, here is what using a Python package in R might look like:
## setup
Sys.setenv("RETICULATE_PYTHON"="/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3")
use_python("/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3", required=TRUE)
reticulate::py_config()
scr <- import(module="scrublet", convert=FALSE)

## usage
E150_scrub <- scr$Scrublet(r_to_py(seurat_object@assays$RNA@counts)$T$tocsc(), expected_doublet_rate=0.06)
doublet_results <- E150_scrub$scrub_doublets(min_counts=2, min_cells=3, min_gene_variability_pctl=85)

The Python package function calls look like scr$Scrublet(). There are some very convenient features like easy conversion of numpy tables to dataframes and the like. If you use RMarkdown, you can even display inline plots from matplotlib. I highly recommend giving it a try, it smoothly bridges the gap between the spheres of development of bioinformatic tools.


text:
Why do I have >10,000 cells in the 10X matrix produced by cellranger?
In my scRNA seq experiment, single-cell libraries were generated using the GemCode Single-Cell Instrument and Single Cell 3′ Library & Gel Bead Kit v2 and Chip Kit (10x Genomics) according to the manufacturer’s protocol. About 8700 cells were added to each channel with a targeted cell recovery estimate of 5000 cells.
However, after processing my raw data with CellRanger software, I get 737,280 cells (barcodes) and about 31,000 genes. I cannot understand why I get so many cells and why this barcode whitelist is not up to 10,000 cells. After filtering this data retaining genes expressed in >= 5 cells as well as retaining cells with at least 200 detected genes, I get around 2200 cells.
Thanks

The raw data from cell ranger contains all of the barcodes detected in the experiment.  These raw barcodes are filtered by cell ranger to identify the barcodes that likely represent cells rather than empty droplets / dead cells. You can find the filtered data in the filtered feature-barcode matrices.  You can also examine the Run summary HTML file to see a summary of the cell ranger barcode filtering.
The reason why there are so many raw barcodes is an important optimization by 10X to increase the cell capture rate compared to techniques such as drop-seq.  In a 10X experiment the number of gel beads is much higher than the number of cells which leads to a higher probability of a cell and a bead ending up in the same droplet.  This also leads to a large number of empty droplets with just a bead and some ambient RNA or a dead cell.  These all end up being sequenced and the empty drops tend to have a small number of unique molecules (the majority less than 10 UMIs / barcode).


text:
How I can reproduce this heat map
I have single cell RNA-seq on 8 time points. I have clustered each time point to 2-3 clusters by seurat. I also have a list of differentially expressed genes between all my 1562 cells in 8 time points. I know how to plot this list of genes in heat map on each cluster of cells or a big heat map of my 1562 cells, but I have also sub-clusters for each cluster (time point).
Could you please help me to figure out how can I have a heat map with sub-clusters for each time point on which I am showing a list of differentially expressed genes?
Something like  the picture below where each day has been clustered to 2-3 sub-population then they plotted common differentially expressed genes between day points
Edit
Lets say 500 genes commonly are being differentially expressed in 8 time point. time point 1 has 3 sub-clusters of cells and time points 2 has 2 sub-clusters of cells. time point 3 and 4 both have 2 sub-clusters of cells. I want to heat map these 500 genes in cells in time point 1 (3 sub-clusters), time point 2 (2 sub-clusters), time point 23 (2 sub-clusters) and time point 4 (2 sub-clusters) in one unified heat map on which cells in each subclusteres (a,b and c) are in column and genes in rows but I don't have any idea how to do that. column annotation function in ComplexHeatmap likely works in this regard But I don't know how to define cells inside each  sub cluster in different colours
I can heat map of 500 genes only on time points easily by seurat but without underlying sub clusters
> DoHeatmap(object = SubsetData(object = object, max.cells.per.ident =
> 100), genes.use = 50 genes, slim.col.label = TRUE, group.label.rot =
> TRUE)

Figure 2 in this article.

Some sample expression data to plot in said heatmap will be of use to encourage answers to your question.

What have you tried to do? Do you have each heatmap for each day and with the same order of genes?

I have independently clustered cells in each time point by seurat. I also obtained differentially expressed genes in each time point separately so I have common genes differentially expressed genes between each pairs of neighbour time points.

In case it is still helpful as the post is rather old, the code below would generate a heatmap with annotations thanks to the ComplexHeatmap package.
But before that I would like to stress out the importance of example data, I bet a lot of experts could not help you just because of they did not have single cell data at their disposal.
I am using some scRNA-seq data stored within a Seurat object. The expression data come from the @scale.data slot and the annotations are from the @meta.data slot. I am using a few general cell type markers as genes and 3 columns from the metadata slot (cell type, cell subtype, sample region). When using the ComplexHeatmap package, the trick is to have the same order in the annotation and actual data tables and this is how you should control what (in your case single cells) is plotted in what order.
library(Seurat)
library(dplyr)
library(ComplexHeatmap)
library(circlize) # for the colorRamp2() function

seurat_object <- readRDS(...)

# Annotations

meta_data <- seurat_object@meta.data

# order of annotations/colors are defined here
ordered_meta_data <- meta_data[order(meta_data$cell_type), ]

# OPTIONAL: YOU CAN PICK COLORS FOR EACH LEVEL OF ANNOTATION
# HERE I PROVIDE COLORS FOR ONLY TWO OF THE THREE LEVELS, COLORS FOR THE
# REMAINING LEVEL IS TAKEN CARE OF BY THE PACKAGE
annotation_colors <- list("cell_type"= c("Epithelial" = "red",
                                         "Immune" = "orange",
                                         "Other" = "gray",
                                         "Stroma" = "green"),
                          "sample_region" = c("B" = "blue",
                                              "C" = "yellow",
                                              "N" = "maroon",
                                              "T" = "darkblue"))
ha = HeatmapAnnotation(df = ordered_meta_data,
                       show_annotation_name = TRUE,
                       col = annotation_colors)

# Expression data

genes_to_use <- c("EPCAM",
                  "CD3D",
                  "NCR1",
                  "CD79A",
                  "LYZ",
                  "COL1A1",
                  "CLDN5",
                  "S100B",
                  "KIT")

seurat_object <- ScaleData(seurat_object,
                           genes.use = genes_to_use)

my_data <- seurat_object@scale.data

# COLUMN ORDER OF THE EXPRESSION DATA SHOULD MATCH THE ROW ORDER OF THE
# ANNOTATION TABLE
my_data <- my_data[, rownames(ordered_meta_data)]

# Heatmap

col_fun = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
png(filename = "example_heatmap.png",
    width = 1000,
    height = 1000)
Heatmap(
  my_data,
  col = col_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  column_order = NULL,
  show_row_dend = FALSE,
  show_column_dend = FALSE,
  show_row_names = TRUE,
  show_column_names = FALSE,
  use_raster = TRUE,
  raster_device = c("png"),
  bottom_annotation = NULL,
  top_annotation = ha
)
dev.off()


text:
Finding genes especific to microenvironment
I have RNA-seq .bam files for 3 patients, tumour and its matched derived model, namely organoid, but I don't have any matched normal sample.
Differentially expressed genes between a tumour and its model mainly come from the presence and lack of microenvironment in cancers and organoids respectively and response to stimuli (artificial environment in organoid culture vs cancer.
How I can get the differentially expressed genes specific to microenvironment?
Any intuition please?

Sauret package for unsupervised learning is a start, lots of discussion on this board about it

I can't find a package named Sauret, if you refer to Seurat it is for scRNA-seq which this question is not tagged with it so I think it won't work well in this case.

The genes related to the environment can be identified using a linear model where each organoid is compared to the correspondent original tumor. This is assuming that the environment for all the organoids is the same and that the differences between tumors is just from different patient but are from the same organ/region.


text:
Discordance of Gene Ontology data between Mouse Genome Informatics and Ensembl sources
In my first approach to the use of Gene Ontology, I have found a discordance if we consult either Mouse Genome Informatics (MGI) or Ensembl for these data.
Concretely, for the mouse genes Zic1 (NCBI 2271) and Zic2 (NCBI 2272), MGI assigns them  up to 23 and 24 different GO terms, respectively.
However these numbers are different according to Ensembl data retrieved through R/Bioconductor. In particular, for Zic2 (Id 22772) the difference is very significant (24 terms in MGI vs. only 3 terms in Ensembl):
> my_ensembl = useMart ("ensembl", dataset = "mmusculus_gene_ensembl")
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22771', mart = my_ensembl)
        go_id
1  GO:0003676
2  GO:0046872
3  GO:0003677
4  GO:0005634
../..
22 GO:0008589
23 GO:0007389
24 GO:0021510
25 GO:0042472
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22772', mart = my_ensembl)
       go_id
1 GO:0003676
2 GO:0005634
3 GO:0016604
>

So, my questions are: Which is the reason for this discordance? Has it something to do with the different types of evidence to assign GO terms? Also, in order to proceed to analysis like GO enrichment, which of these two databases would be more appropriate?

I'm looking into this. We should have all the annotations from MGI and I'm not sure why we don't.

Short answer, we get our annotations from GOA, which has the three listed. We're in touch with GOA to see why they differ from MGI.


text:
Error Code in WGCNA: Dealing with large data sets: block-wise network construction andmodule detection
I am comparing single block module detection with 2 block module detection in WGCNA.
Does anyone know what this relates to? I have rechecked and rechecked my codes for both the single and 2 module detection and cannot see an error. Thank you
I get the error code:
Error in plotOrderedColors(dendro$order, colors = colors, rowLabels = rowLabels,  :
  ERROR: length of colors vector not compatible with number of objects in 'order'.

Using the following code:
net= blockwiseModules(Th1_h5, power=6,
                      TOMType = "signed", minModuleSize=30,
                      reassignThreshold =0, mergeCutHeight = 0.25,
                      numericLabels = TRUE, pamRespectsDendro = FALSE,
                      saveTOMs = TRUE,
                      saveTOMFileBase = "PBCreponseTOMB",
                      verbose= 3
mergedColors = labels2colors(net$colors)

plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

bwet = blockwiseModules(Th1_h5, maxBlockSize = 8000,
                        power=6, TOMType = "signed", minModuleSize = 30,
                        reassignThreshold = 0, mergeCutHeight = 0.25,
                        numericLabels = TRUE,
                        saveTOMs = TRUE,
                        saveTOMFileBase = "Th1_NR_R_TOM-blockwise",
                        verbose=3)
bwLabels = matchLabels(bwet$colors,moduleLabels)

bwModuleColors = labels2colors(bwLabels)

plotDendroAndColors(geneTree,
                    cbind(moduleColors, bwModuleColors),
                    c("Single block", "2 blocks"),
                    main = "Single block gene dendrogram and module colours",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang =  0.05)

You will need to provide more detail. How big was the data from which you calculated the network (how many genes) and how did you determine geneTree?

I have 12711 genes after filtering and have added the additional code above. It has worked with the simple Automatic, one-step network construction and module detection and Step-by-step network construction and module detection but not when trying with 2 blocks. Thank you kindly

The problem is that blockwiseModules split your data into (probably 3) blocks, because the default maxBlockSize is 5000 and smaller than the number of genes in your data. If you have enough RAM (16GB is recommended for 12k genes, but 8GB just might be enough), rerun blockswiseModules again with maxBlockSize raised to above the number of your genes. Since you were able to run the function with maxBlockSize=8000, hoepfully you'll be able to do the single block calculation.

Thank you so much! I will give it a try now, thank you!

It works! Thank you so much, as always for your time and help!

authors: Кристина Павленова


text:
How can I convert Affymetrix microarray data to bedgraph format?
I've previously asked a different question about the ChIP-chip microarray data that I'm working on. I'm new to this so I'm sorry for asking here again.
I have .CEL files and I'm using the Affy package for R to read them. After I load all the files in, I then normalise the data and  generate a .txt file containing the expression values of each file (ranging from 11 to 12). Column 1 contains the gene names and the rest of the  columns contain the expression values for each of the files.
I was wondering if there is a way to convert these data to bedgraph or BED format? What I'd like to do is a metagene analysis but I don't know how that would be possible with expression value data.
I looked online and there are some kinda relevant suggestions/solutions, but they are for Windows machines and I'm working on MacOS.

Do the generated txt files contain the physical coordinates of the probes? If so, converting them to BED should be a simple matter of rearranging the columns (or they might already be valid BED files).

Hi Konrad, thank you for your response.
The generated txt file does not contain the coordinates of the probes, however, I have a gff file that contains all this information and I guess I could look up the coordinates for each gene of the array. However, I don't know how the expression values could be converted?

Why makes the solutions for Windows specific for Windows? Did you try to overcome these using equivalent programs or formats?


text:
Chunk alignment in a name sorted bam for parallel processing
I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.

Could you confirm what you mean by "chunks". Do you mean a genetic locus, e.g. X nucleotides is one "chunk". Alternatively, do you mean subsets of the read depth?

I want to split the alignment into small segments, e.g., 10M alignments per segment. The expected input bam file is name sorted, so roughly 1-10M alignments in the first chunk,10M - 20M alignments in the second chunk and so on. The problem is that I don't want to put the alignments from the same read into different chunks. Thank you.

What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing.

What have you tried? In particular, what python code using pysam have you written to do this? The steps are mostly (A) name sort and (B) iterate over and write a lot of files.

I think it is a bit outside of the scope of this site to write a program for you. Perhaps if you show us what you have tried so far, it would be more possible for us to assist you. In the question's current state, you are asking us to write you a program which incorporates pysam.

Thank you all for your help. Did some search. I think that for fastq files or sam files, it should be easy based on the class here: https://www.blopig.com/blog/2016/08/processing-large-files-using-python-part-duex/
However, for bam files it's more challenging as I don't know how to use 'seek' in pysam. Or maybe just treat bam files as special zipped files and use general gzip reader for this purpose? Or is it possible to do so by using pysam? Thank you all.

Basically, I want to know where all the alignments of a read mapped to, so I decide to process read by read. Here I distinguish reads from alignments.

Which aligner was used? bwa stored alternative mapping positions in the XA tag.

These are RNA-seq data and I use STAR for alignment. I can see four tags, e.g., NH:i:6 HI:i:6 AS:i:49 nM:i:0

You can sort by read name instead of position. This is relevant for whatever you are doing. You should read this thread to understand the bit flags for multiple primary alignments. It should be helpful.

I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing code if going the be the slow bit.
import pysam
infile = pysam.AlignmentFile("input.bam")

chunk_size = 10000000
outfile_pattern = "output_segement%d.bam"

chunk = 0
reads_in_this_chunk = 0
old_name = None
outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

for read in infile.fetch(until_eof=True):

    if old_name != read.query_name and reads_in_this_chunk > chunk_size:
        reads_in_this_chunk = 0
        chunk += 1
        outfile.close()
        outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

    outfile.write(read)
    old_name = read.query_name
    reads_in_this_chunk += 1

outfile.close()

Thank you, this works great! What I can think of to optimize the code a little bit is 1) I don't need to check the name of each alignment. Only after reading a specified number of alignments, I check the read name, read again until I find a read with a different name. 2) It seems un-doable for pysam to pass the start and stop positions of chunks to multiprocessing.pool.

Not sure how can I accept this answer.

@jenny00 the comparison of read names will take next to no time. The slow bit is reading the read into python, once its in manipulations are fast.
And as for 2: No, pysam is not compatible with python's simple multiprocessing API because you can only pass picklable objects through to the child processes and AlignedSegments (or reads) are binary objects - they are actually a thin adaptor layer on the C++ based htslib object.

Thanks, I just noticed that. So I have to put the necessary information from each alignment into a dictionary and pass dictionaries down for parallel processing. It does not save much time compared to just using one process.

If you are using pythons multiprocessing.Pool infrastructure, you might find that the overhead associated is high - each submission to a Pool causes the creation of a new python process, which does its processing and then the process is closed down (i think). You might like to consider the more complex to implement asyncronous worker and processing queue type set up. Also try increasing the size of the chunks - the fewer chunks, the small the number of process spin ups and close downs.

Thanks. Increasing the chunk size and using map_async help. But it is only twice faster using four cores, compared to using a single core.

python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
  BEGIN { part=0; line=n }
  /^@/ {header = header$0"\n"; next;}
  { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" }
    else { print $0 >> part".sam"; line++; }
    last_read = $1;
  }'

The output will be in sam format.
You can controll the number of reads per chunk by setting the n value for the awk command.
EDIT:
A second version that will output in bam:
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
BEGIN { part=0; line=n }
/^@/ {header = header$0"\n"; next;}
{ if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}
  print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam"
  last_read = $1;
  line++;
}

Thank you so much! This works great as well, but I want to use pysam so the code can be integrated into a longer function. I don't even need to save the chunked alignments, but just pass them for parallel processing. It could be quite useful if I want to change my design to split and save a bam file into smaller bam files. I think that many other people may find this solution helpful as well.

I was considering awk too. The problem with Perl/Python is the code will process the sequence data "[horizontal] line by [horizontal] line", when awk can do this "[vertical] line by [vertical] line "


text:
Total reads aligning to each reference within a bam file
I have two PCR amplicons that have been multiplexed and sequenced using the nanopore minion.
I have aligned the fastq reads using minimap2 with a reference file containing both amplicon sequences and generated a bam file that I have viewed using IGV.
I am looking for a way to generate some simple summary statistics.
In particular, is there a way to extract the total number of fastq reads aligning to each amplicon reference from the bam file?

The quick way to get the number of alignments on each reference is
samtools idxstats my_bam.bam

Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.
In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:
samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam

This one-liner below will work better for long reads than samtools flagstat in that it only counts the primary alignment for each read and samtools flagstat doesn't seem to calculate some stats for long reads. I have never seen samtools flagstat output stats on a per-reference basis, but am curious if so!
This answer filters out secondary and supplementary alignments for your reads (-F 2304) that have some alignment to both amplicon reference and just keep the best one. This might give a more accurate idea of how many reads of each amplicon are in the library.
samtools view -F 2304 myfile.bam | awk -F $'\t' '{a[$1, $3]++} END{for (i in a) {split (i, sep, SUBSEP); print sep[1], sep[2], a[i]}}' | uniq | awk '{print($2)}' | uniq -c | sort -k1 -nr

adapted from this.

This can be done using
   samtools flagstat your_bam_file

I think that flagstat counts the number of alignments in the bam file.


text:
Mapping RefSeq-ID or NCBI-ID (Entrez-ID) to EC number
I annotated my bacterial genomes using the new NCBI Prokaryotic Genome Annotation Pipeline and now, I want to annotate EC-numbers.
In the master annotation file (.gff), I get the RefSeq-ID of every successful annotation, e.g. RefSeq:YP_805528.1.
When I visit the respective NCBI page, I see that it sometimes is mapped to an EC number.
I'd like to automate this process: Convert every RefSeq-ID to an EC number. How can I do this?
(I may already have discovered the first step: I can map the RefSeq-ID to the NCBI-ID (Entrez-ID) using this NCBI conversion database: YP_805528.1 -> 4419225. But I don't see how to get from these databases to EC.)
Are there other ways? As far as I know, most people blast their genes to already-mapped databases like SwissProt and ExPASy ENZYME. But I don't like this approach because technically, the NCBI annotation pipeline already did this. Why compute the same thing twice?

You can use Entrez Direct for this as follows:
efetch -db protein -id YP_805528.1 -format gpc -mode xml | xtract -insd Protein EC_number
YP_805528.1     6.3.2.13


text:
Multiple correlation with R
I have a data frame presenting different variables in column and for different species name in rows as follow. In the column a have only one independent variables(Proteome size) and all the others are independent variables. I want to use the PIC package to correlate all these dependent variables separately to my independent variables, a way to generate the R coefficient and p-value.
data <- read.table("t1.txt", row.names = 1, header = T,sep='')

data
path1 path2 path3 path4 Proteome_size
ahli 4.0 -3.2 0.4 -1.0 -0.2
alayoni 3.8 3.4 -1.8 2.2 -0.3
alfaroi 3.5 2.6 0.8 -2.4 -0.5
aliniger 4.0 0.1 -1.7 2.4 0.0
allisoni 4.4 2.0 -3.7 0.5 -0.1
allogus 4.0 -2.8 0.6 -1.0 0.6
altitudinalis 3.8 2.9 -6.1 2.3 1.2
alumina 3.6 0.7 1.5 -2.7 -1.3
alutaceus 3.6 1.2 -0.6 -1.7 0.3
angusticeps 3.8 4.6 -2.0 1.2 1.7

I have this code following code that can correlate any one of the ‘path*’ column to proteome size:
data <- read.csv("data.csv", row.names = 1
tree <- read.tree("anolis.phy")
Proteome <- data[, "Proteome_size"]
awe <- data[, "path1"]
names(Proteome_size) <- names(var1) <- rownames(data)
hPic <- pic(proteome_size, tree)
aPic <- pic(gene1, tree)
picModel <- lm(hPic ~ aPic - 1)
summary(picModel)

My question is that, how can I implement this code to do that for all the columns simultaneously and store the statistics info as (R and p_value) for each variable somewhere ?
Thanks

Welcome to the site. Have you tried to do a loop for each column to correlate them with the Proteome_size?

Thanks for your quick response. I am a beginner  in R.  To do that  need to create a vector before  as : gene.names <- colnames(data[, "proteome_size"]). unfortunately it failed and the  vector is NULL with length 0.  I am not so good in for loop too

Sorry, I don't actually understand what the variable awe is doing. Also I don't see which variable will change through the loop. It's quite easy to write a for loop in R if you can tell which variable needs to be investigated!

Your code gene.names <- colnames(data[, "proteome_size"]) produces a NULL vector since R is case sensitive and your data frame does not have a "proteome_size" column.

One option is to use one of the apply family functions. If you enclose your code above (which generates errors by the way) in a function, you can then "apply" this function on the rows or columns of your data with something like apply(X, MARGIN, FUN, …), the margin being 1 for rows or 2 for columns (as in your case).

Like haci mentioned, you can use the apply family of functions to loop over all your columns.
apply(data,2, function(path_name){ summary(lm(data$Proteome_size ~ path_name -1))$adj.r.squared })

this gives you the adjusted r square value for your linear model, and
apply(data,2, function(path_name){ summary(lm(data$Proteome_size ~ path_name -1))$coefficients[,4] })

would give you the p-values.
What happens here is apply takes your matrix and chops it up into pieces(a piece is a column here because I gave the option 2) and each pieces is fed into the function as 'path_name' and the output for all the pieces is given back to you.
You might get a warning about getting a perfect fit because I am correlating 'Proteome_size' against itself.
Also,
awe <- data[, "path1"]
awe <- data$path1

are the same when dealing with columns and the second can be more convenient because you can search for autocomplete by hitting tab after typing the $


text:
Error compiling HMMER v3.2.1
I tried to compile HMMER v3.2.1 in Cygwin using
./configure

followed by
make

or
make install

but failed. I got this config.log file to troubleshoot and I suspect the solution is in this line:
configure:2908: error: no acceptable C compiler found in $PATH

What does it mean? Is it really all I need to address? I followed the installation steps in this document.

It means you need to download a C compiler. I would recommend using gcc!

Have a look here on how to install gcc https://askubuntu.com/questions/271388/how-to-install-gcc-4-8

gcc is a C compiler

I would strongly recommend to use Bioconda instead of compiling. It is one of the channels ("distributors of packages") in Anaconda Cloud which is basically a package source for conda.
In short - conda allows you to create virtual environments isolated from system and to install various precompiled software packages with one command (so, by the way, you can easily use multiple software versions on one system). Also, Bioconda is just one channel aiming on various bioinformatics software, but many more packages can be found in aforementioned Anaconda Cloud.
Tutorial how to start using Bioconda is on its homepage, and it is really easy (far more easier than to struggle with compiling).
Bioconda HMMER package
P.S. Of course there are cases in which you need a specific build of given software. Than it is really needed to compile manually.

Thank you. The website says "Bioconda supports only 64-bit Linux and Mac OSX." I am a Windows user, so I assume I will need to go back to compiling. Bioconda sounds amazing though; thanks for bringing it to my/our attention.

I tried to set up channels with conda config --add channels defaults, conda config --add channels bioconda & conda config --add channels conda-forge in my Anaconda prompt which seemingly worked, but conda install hmmer didn't.

@JernejT On Windows? And what does conda install hmmer output? Anyway, Linux is much more supported platform for bioinformatics. Worth to try :)


text:
How do I interpret my phylogenetic tree construction from Glucosidases in the Actinobacteria class?
Dear knowledgeable bioinformaticians, I'm currently pursuing an undergraduate degree and as a result I do not have so much knowledge in this topic. I'm having some troubles understanding my phylogenetic three construction.
For the construction of a phylogenetic tree UGENE was used. 11 homologous proteins from different species in the same genus/phyla were aligned. For the tree construction PHYLYP Neighbour Joining (NJ) with the distance matrix model “Jones-Taylor-Thornton” was used. Two trees were constructed, one without bootstrapping and one with bootstrapping (1000 pseudoreplicates).
How do I interpret my values? Are they unusable since I have low bootstrap-values in a lot of my nodes? All guidance is appreciated, thank you.

Welcome to the site. Could you please clarify which one is without bootstrapping? It seems to be the one on the top, but I might be wrong

Thank you! The upper one is not bootstrapped, the lower one is the bootstrapped one.

Bootstrap values above 750 in the lower tree represent robust clades. Bootstraps above 800 are preferred and anything above 900 is strongly supported.

Bootstrap values usually give you useful information, be they low or high. Usually, a "high" bootstrap value is something >= 70%. That means a clade is strongly supported; lower values are telling you the clade is not clearly separated from others. This could be for several reasons, e.g.:
1) A subpopulation is not really separated from others, e.g. the viral sequences in  town "A" and town "B" are completely admixed, meaning the population is panmictic and there is continuous gene flow between the two of them;
2) the populations diverged not long ago, it's too soon to have any phylogenetic signal, at least using the marker you are using (different genes evolve at different rates). Usually, housekeeping genes (which are usually very conserved) are used for populations that supposedly diverged long time ago.
The interpretation of your tree should take into account the biology and epidemiology of your target organism. If you really want to learn more, I suggest a lot of reading, and definitely a class dedicated to molecular phylogenetics.
This might be a good start:
https://www.amazon.com/Phylogenetic-Handbook-Practical-Approach-Hypothesis/dp/0521730716
Also, the NJ algorithm should allow for uneven branch lengths, maybe the tree you are showing is the consensus tree?
Hope it helps
Max

The branch lengths are defaulted, however, they do allow different lengths based on how one would like them to look like. It's just a display function in UGENE. Thank you a lot for your comment, extremely helpful!

You're welcome, and good luck. If you found the comment useful, please consider upvoting it or accepting it as answer.

"the epidemiology" infers the bacteria are pathogenic.


text:
Error in as.vector(x) : no method for coercing this S4 class to a vector
I tried to run the following code in R studio. Everything worked fine, except at the last step [write.table(mdat, "recount_mdat.csv")] when I tried to export the 'mdat', I got the following error:
Error in as.vector(x) : no method for coercing this S4 class to a vector

My code:
 library('dplyr')
 library('recount')
 library('magrittr')
 library('ffpe')
 library('RSkittleBrewer')
 library('SummarizedExperiment')
 library('devtools')
 trop <- RSkittleBrewer::RSkittleBrewer('tropical')
 colon_proj <- c('SRP025982', 'SRP42161')
 if(any(!file.exists(file.path(colon_proj, 'rse_gene.Rdata')))) {
 sapply(colon_proj, download_study)
 }

 blood_proj <- c('SRP066834', 'SRP041736', 'SRP060416')
 if(any(!file.exists(file.path(blood_proj, 'rse_gene.Rdata')))) {
 sapply(blood_proj, download_study)
 }
 proj <- c(colon_proj,blood_proj)

 dat <- lapply(proj, function(x) {
 load(file.path(x, 'rse_gene.Rdata'))
 return(rse_gene)
 })
 proj

 sapply(dat, dim)

 metadata <- all_metadata('sra')

 write.table(metadata, "recount_metadata.csv")

 mdat <- do.call(cbind, dat)

 write.table(mdat, "recount_mdat.csv")

is(mdat)
## "RangedSummarizedExperiment" "SummarizedExperiment" "Vector" "Annotated"

What is the class of mdat? Did you try is(mdat)? What did you get?

@IIrs Hi I got this after this step
           [1] "RangedSummarizedExperiment" "SummarizedExperiment"
           "Vector"
   [4] "Annotated"

Great, see how fast you got your answer when you provided enough details, remember to include all the relevant details next time! Glad you found how to solve the problem

You can't cbind a bunch of obscure object types.
If you want merged count tables you should do this:
mdat <- do.call(cbind,lapply(dat,assay))

Where row.names are Ensembl gene IDs and col.names are the SRR accessions.
Then run your table writing command.
-
If you want the coordinates of your genes then do this to make a bed with the genomic locations:
library(rtracklayer)
export.bed(rowRanges(dat[[1]]),'gene_coords.bed')


text:
Reconstruct a protein's 3D structure from a pair-wise distance matrix
What is a good method to reconstruct a protein's 3D structure from a pair-wise distance matrix of residues?
I think the question is self-explanatory. I know of this method
but it seems that it uses contact maps, which actually use a certain cut-off. I don't have a cutoff based contact matrix but a matrix with actual distances.

You have distance between alpha carbons or beta carbons (or other atoms of the residues)? Usually a threshold is set because below certain measure it can't be and above certain threshold it doesn't mean a contact between the residues. Do you have any other information besides the distance matrix, like the sequence of the protein?

@llrs They are only carbon alpha distances and I want to reconstruct the backbone with carbon alphas and no more atoms.No more information, just the distances between carbon alphas for each residue.

You could use the multidimensional scaling approach, but I suspect it might propose invalid angles between the backbone. It would be possible to use it as a starting point before checking the angles between the backbone. But I don't know if there is any software to do this (it is a strange problem to have), or either you know the sequence or the sequence and the structure, but just the distance is strange.

Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.
But I don't enough about the context of your problem to tell if/how it is applicable.


text:
Calculating alignment/mapping time
I  am trying to assemble a plant genome using AWS resources using velvet. Plant genome is huge (> 10 times human genome) and coverage is around 30 x. We are planning for de novo assembly with Velvet (threads enabled).
I would like to know if there is a calculator that can provide approximate time taken for assembly, given the necessary details. For eg. if I furnish RAM, CPU (type of instance), genome size, approximate coverage, type of sequencing (PE or SE) and number of client nodes, it should give me the number of hours or days that would take for assembly (reference and/or de novo).
Note: Previously posted on biostars, but no replies.

I don't know much about assemblies, but do you know how complex it is? chromosomes, repeated regions and other information about the complexity of the DNA? Are you using long reads or short reads (I think that short-reads, but just to make sure). Lastly, could you provide the link to the post on biostars? (Just so if someone answers there people who find this post know where to look for)

BioStars link: https://www.biostars.org/p/348943/

IMO this is unknowable, which is why you received no reply on biostars.

Velvet won't work for a 30Gb genome for sure. You might have a better chance with assemblers capable of human assembly, such as soapdenovo and discovar. But even with them, probably you won't get reasonable results.