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:
Why does repeatmasker annotate transposons of length 1?
I am working on Mac OSX. I have downloaded RepeatMasker from www.repeatmasker.org. I have downloaded RepBase from www.girinst.org. I have download the WS266 version of the C. elegans genome from WormBase. I have installed nhmmer using conda.
I have then run the following command:
repeatmasker -s -lib RepBase23.08.fasta/celrep.ref -dir annotation/
-gff c_elegans.PRJNA13758.WS266.genomic.fa

I then imported the annotation table into R. I notice that some transposons have a length of 1!
How is this possible?!
Also, how can repeatmasker calculate such percent divergence, deletion, insertion with a query width of 1? For some of these predictions, both the query and repeat lengths are 1 (so an alignment of 1 nucleotide against 1 nucleotide!)! But the Smith-Waterman scores don't match this. There is something important I am missing here. Can someone help me interpret the table below and explain why some transposons have a length of 1 in the genome?
Thanks.

The image only contains transposons with a width of 1 as far as I understand. I've selected the rows where queryBegin == queryEnd. The queryBegin variable = queryEnd for all the lines here. It is a subset of the original data frame.

source: Stack Exchange

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

sort: votes


text:
Using ONT MinION, why is that we cannot get a full length DNA read?
Here, https://www.youtube.com/watch?v=E9-Rm5AoZGw, ONT claims they can read non-fragmented DNA and that with the hairpin like structure the double strands can be read at once. Then why is that MinION reads are of variable in length and not the length of a whole DNA? BUt here:https://www.youtube.com/watch?v=e9wtAIvPPxY, ONT claims that the read lenght can be configured, what does this mean?

The MinION will read the full length of the the fragments which pass through the pore. Are you asking why it can't read the full genome in one go?

The MinION should be able to sequence DNA fragments entirely, regardless of the size, if you can get them to the pore intact. That is the problem: during DNA extraction and library preparation you will introduce breakage.

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


text:
Visual representation of detecting deletion by different read types
I'm having some trouble differentiating between CNV, split reads and read pairs and how each one elucidates the presence of a structural variation, say deletion.
Given a reference and donor sequence, I'd like to see a visual depiction (a link to a reference would be welcome) of a deletion event being detected by each of the above read types. In particular, I'd like to see how each of the read types detects the deletion event itself (i.e structural variant type) and its location (i.e the breakpoint junction).
Bonus points if it also includes some of the other common structural variants like inversion, duplication and insertion.

Welcome to Bioinformatics! This question need more information to be possible to answer. What papers did you read or which tools aren't helpful (and why not?)? Many thanks!


text:
What does that mean if there is no bifurcations in my time series single cell data?
Sorry,
I have analysed my time series single cell data with Scuba algorithm but that says that there is no bifurcations in my data and outputted a single chain of 4 states with no branching like below figure

Sorry, I googled but I am not quite sure what this means
Any suggestion please?
Thanks a lot

The most obvious interpretation is that the four cell lines are related in a linear pathway rather than a branching one.


text:
Group genes by functional categories suming expression values
Using the count of rpkm values from genes in a metagenome sample, I want to group these genes into established categories (for example KEGG or COG). For each sample, my goal is to determine which categories are better represented in each one.
Considering the example table whit genes, samples (1-5) and categories.
Each gene has their expression and the column categories classify them.
import sys
import numpy

dictionary = {}

table = open(sys.argv[1], "r")

for line in table:
    next()
    cols = line.strip().split("\t")
    gene = cols[0]
    s1 = float(cols[1])
    s2 = float(cols[2])
    s3 = float(cols[3])
    s4 = float(cols[4])
    s5 = float(cols[5])
    counts = [s1, s2, s3, s4, s5]
    kegg = cols[6].strip().split(" ")[0].replace('"','')

    if kegg not in dictionary:
        dictionary[kegg] = [counts]
    else:
        dictionary[kegg].append(counts)

#print dictionary
for k, v in sorted(dictionary.iteritems()):
    m = numpy.array(v)
    print k,  [sum(m[:,i]) for i in range(5)] # range(number_of_samples)

Original table to count the columns values:
gene    s1  s2  s3  s4  s5  Category
name01  0   2   2   0   0   A
name02  3   0   1   0   0   A
name03  0   2   1   0   0   A
name04  0   0   1   0   0   B
name05  5   0   1   0   0   C
name06  1   0   0   0   0   D
name07  2   0   0   0   0   D
name08  0   0   3   0   0   E
name09  1   0   0   0   0   F
name10  1   0   0   0   0   F
name11  3   0   0   0   0   F

The result obtained is this:
"""
# eliminated from the question

Cat s1  s2  s3  s4  s5
A   0.0 2.0 1.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   2.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   3.0 0.0 0.0 0.0 0.0
"""

But my expected (obtained) table is this:
Cat s1  s2  s3  s4  s5
A   3.0 4.0 4.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   3.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   5.0 0.0 0.0 0.0 0.0

Have you looked into using the pandas library?  It is specifically designed to do these transformations.

not yet. Thank you for your suggestion. If you have some idea how to use it, I thank you

I suspect you do not build your dictionary correctly. If you append a list to a list, you get a list inside a list. extend is probably what you want. See https://stackoverflow.com/a/252711/1878788 Besides, you might be interested in collections.defaultdict: dictionary = defaultdict(list) and then: dictionary.extend(counts), without need to test whether it already has a kegg entry.

@bli I corrected the script and updated the question. Now it works.

As suggested in the comments, the pandas module is quite convenient for this type of work.
Here is how you could do (assuming you have tab-separated values as input, and you want similar file format for the output):
import pandas as pd

counts = pd.read_table("counts.tsv")
# At this point, counts looks as follows:
#       gene  s1  s2  s3  s4  s5 Category
# 0   name01   0   2   2   0   0        A
# 1   name02   3   0   1   0   0        A
# 2   name03   0   2   1   0   0        A
# 3   name04   0   0   1   0   0        B
# 4   name05   5   0   1   0   0        C
# 5   name06   1   0   0   0   0        D
# 6   name07   2   0   0   0   0        D
# 7   name08   0   0   3   0   0        E
# 8   name09   1   0   0   0   0        F
# 9   name10   1   0   0   0   0        F
# 10  name11   3   0   0   0   0        F

summed_by_cat = counts.groupby("Category").sum()
# At this point, summed_by_cats looks as follows:
#           s1  s2  s3  s4  s5
# Category
# A          3   4   4   0   0
# B          0   0   1   0   0
# C          5   0   1   0   0
# D          3   0   0   0   0
# E          0   0   3   0   0
# F          5   0   0   0   0

# Change the index name to get the desired first column header in the output:
summed_by_cat.index.name = "Cat"

# Write to a file
summed_by_cat.to_csv("by_cat.tsv", sep="\t")

And resulting the file contains the following:
Cat s1  s2  s3  s4  s5
A   3   4   4   0   0
B   0   0   1   0   0
C   5   0   1   0   0
D   3   0   0   0   0
E   0   0   3   0   0
F   5   0   0   0   0

You could also load your table using counts = pd.read_table("data.tsv", index_col="gene") and counts would have the gene names as index instead of default numbers:
        s1  s2  s3  s4  s5 Category
gene
name01   0   2   2   0   0        A
name02   3   0   1   0   0        A
name03   0   2   1   0   0        A
name04   0   0   1   0   0        B
name05   5   0   1   0   0        C
name06   1   0   0   0   0        D
name07   2   0   0   0   0        D
name08   0   0   3   0   0        E
name09   1   0   0   0   0        F
name10   1   0   0   0   0        F
name11   3   0   0   0   0        F

This doesn't seem to affect the final outcome.

The complete answer by far.

Below is code to do that, though if these are really gene counts then I highly suggest you provide more details about why you're trying to do this before you actually do it.
#!/usr/bin/env python
print("Cat s1  s2  s3  s4  s5")
v = []
category = None
for line in open("counts.txt"):
    if line.startswith("gene"):
        continue
    cols = line.split()
    if cols[-1] == category:
        for idx in range(5):
            v[idx] += int(cols[idx + 1])
    else:
        if category:
            print("{}   {}".format(category, " ".join(["{:.1f}".format(x) for x in v])))
        v = [int(x) for x in cols[1:-1]]
        category = cols[-1]
print("{}   {}".format(category, " ".join(["{:.1f}".format(x) for x in v])))


text:
Error correction within the long read
I am going to get some data from plasmid sequencing to identify SNPs on the plasmids. What is done in the lab is the following:

The plasmids are purified by size.
We amplify the plasmids using the phi29 polymerase. The polymerase will go through the plasmid multiple times. Hence we get the same sequence concatenated multiple times. My question is related to this step:
Finally, We sequence it using Oxford Nanopore.

My question is:
On step two, I wrote that "we get the same sequence concatenated multiple times". For me, this a potential source of information to correct the base calls prior to aligning them to the reference. Since you have the same sequence multiple times (concatemers).
What I would like to know is if there is a tool that uses the information of the concatemers to improve the base calls. I have tried to find some methods but found none. I know that PacBio has a similar flavour using the "circular consensus calling", but I have not found any methodological explanation.
thanks!

Nice project! Do you know the length of the (linearized) plasmids? I haven't heard of any tool similar to this but I don't usually analyze long-reads sequencing data. Good luck with your project

I don't, the project aim is to detect which plasmids a bacteria contains. Thanks for your answer

Then the project is more complicated, because if you don't know the length of the repetitions of the plasmid you'll need to search for a range of size of the repetitions, which can be longer/more difficult (depending on how you compare the different portions of the reads)

Yes, there is a tool to do this called R2C2 by the Vollmers lab at UC Santa Cruz.
https://www.biorxiv.org/content/early/2018/06/04/338020

That was definetively the tool I was looking for. Thank you very much!


text:
Retrieve all genes for a specific gene ontology
I'm having a bit of confusion about how can I retrieve all the genes that are associated to a specific gene ontology.
Lets say I need to find all the genes that are associated with this GO entry:

GO:0030098 lymphocyte differentiation

This is the graph representing the parents of my GO. In the same page I find the Child Terms that represent the children of my term.
Now, I can download the goa file (related to the organism I'm interested into) containing the GO associations and grep my term to extract a list of associated genes but, as far as I can understand, this list does not represent the whole graph that actually describe my GO term.
Should I also get the genes in the childred terms of my term? Or in the parents terms? If yes, then again I have to go back in the graph and find the parents/children of the found parents/children and so on, till the whole graph is visited.
Is this the correct approach? Or should I only focus on those genes related to my GO:0030098, and that's it?
I found this:
Retrieving a list of human genes having GO associations
and this but I think they don't address my question since they do not describe the parent-children association.
Also I found this, on the GO FAQ, but what I need to do is quite the opposite.
Summarising: I need to find the genes associated to a specific GO term. Should I also look for its parents or children or not?
Any clues?

You can use BioMart to filter by your GO term and get the genes as attributes. BioMart is ontology-aware so will pull out all genes associated with your term and with any of its child terms.
There is no need to look up the GO child terms as BioMart already deals with this. A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.
If you're not already familiar with BioMart, there's a help video.
Here's a BioMart query that gets you exactly what you need in human using the web interface
If you want to work with the R interface to biomaRt, here's some documentation and this is the identical query to above:
library(biomaRt)
ensembl=useEnsembl(biomart="ensembl", dataset='hsapiens_gene_ensembl')
getBM(attributes = c('ensembl_gene_id','external_gene_name', 'go_parent_term', 'name_1006'),
      filters = 'go',
      values = 'GO:0030098',
      mart = ensembl)

ok, so as you're suggesting what I should look, in GO, is for the children of my GO term. But a child of my term can contain genes from other parents, isn't it?
according to this link http://www.geneontology.org/page/ontology-structure it seems that some kind of association may appear between different branches.

No, don't worry about looking in GO. Just put the one term you have into BioMart. It will find all genes associated with your term and all genes associated with its children.
A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.

Here's a BioMart query that gets you exactly what you need in human

@Emily_Ensembl please [edit] the comments to your answer (the comments might be deleted and there are very useful comments). Thanks for answering this

Why would the comments be deleted?

Emily, one more detail. The example query you sent, is equivalent to this piece of code (changing the mart and the GO entry)?
http://127.0.0.1:19280/library/biomaRt/doc/biomaRt.html#retrieve-all-entrezgene-identifiers-and-hugo-gene-symbols-of-genes-which-have-a-map-kinase-activity-go-term-associated-with-it.

I don't have permission to connect to that. Can you send it as a proper URL.

@Emily_Ensembl Because in Stackexchange the are for clarifications and improvements of questions and answers, and as such the content is assumed to be added to the relevant question or answer. Also they are considered "second class citizens"...

Emily, if you click on it then it does not work, if you copy-paste in a new tab, definitely it will. Otherwise, from the page that it opens, I'm focusing on the section 4.6  Retrieve all entrezgene identifiers and HUGO gene symbols of genes which have a “MAP kinase activity” GO term associated with it.

@gabrielet the first section of your link (http://127.0.0.1:19280) is for a local port of your machine! We can't access your machine, so we cannot see the code you are mentioning

ops! sorry Llopis and Emily, you're right. This should work
https://bioconductor.org/packages/release/bioc/vignettes/biomaRt/inst/doc/biomaRt.html#retrieve-all-entrezgene-identifiers-and-hugo-gene-symbols-of-genes-which-have-a-map-kinase-activity-go-term-associated-with-it.

Yes, that's it. But change the attributes to whatever you want to see in your results table.


text:
How to convert data in gmt format to dataframe?
I downloaded c5: gene ontology gene sets file from http://software.broadinstitute.org/gsea/downloads.jsp
I opened the "c5.all.v6.2.symbols.gmt" file in csv format and It looks like below:

I want to convert the .gmt file into a dataframe. And it should look like below:

can anyone say how I can do that. Thanq

This might be an XY problem (you are asking about a problem, but you really have another one). Why do you need this data in data.frame format? Also consider that the GO data is a directed acyclic graph, which has a special structure that should be taken into account

for some analysis I want the data to be looked like the output I posted.

Also note that the GO annotations imported by MSigDB are quite old. If you are interested in associations between GO and genes (updated) you should use other sources

other resources like?

Like the GO database directly or PANTHER. Or to retrieve the results biomart or its package biomaRt

Ok. Answering my own question.
I did this way. And I got the output I want.
install.packages("msigdbr")
library(msigdbr)

m_df = msigdbr(species = "Homo sapiens", category = "C5")
head(m_df)

Output:
gs_name gs_id gs_cat gs_subcat human_gene_symb~ species_name entrez_gene
  <chr>   <chr> <chr>  <chr>     <chr>            <chr>              <int>
1 GO_14_~ M184~ C5     MF        AANAT            Homo sapiens          15
2 GO_14_~ M184~ C5     MF        AKT1             Homo sapiens         207
3 GO_14_~ M184~ C5     MF        ARRB2            Homo sapiens         409
4 GO_14_~ M184~ C5     MF        BAD              Homo sapiens         572
5 GO_14_~ M184~ C5     MF        DAB2IP           Homo sapiens      153090
6 GO_14_~ M184~ C5     MF        DDIT4            Homo sapiens       54541
# ... with 2 more variables: gene_symbol <chr>, sources <chr>

found this solution from here https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html

perfect! thanks.

You can import the csv file as a table. This does not require further libraries:
#you're reading a csv file, using tab as field separator and considering the first line as the headers of the data table
data <- read.csv(filename, sep="\t", header=TRUE)

then what you need are just two columns of your data:
colnames(data) #will give you the names of the columns that can be accessed
#to access two columns and merge them into a selectedColumns table
selectedColumns <- cbind(data$colName1, data$colName2)

The question is tagged with R so it is an R data.frame.

you're right! I missed it

You can use GSA.read.gmt function from GSA package. The following code can be used to convert the file to a dataframe. Just ignore the warnings.
Original_response
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
gene_names <- unlist(data$genesets, use.names=FALSE)
your_dataframe <- cbind(data$geneset.names,gene_names)
colnames(your_dataframe) <- c("Pathways","Genes")
mydataframe <- as.data.frame(your_dataframe)
head(mydataframe)
                                               Pathway   Gene
1            GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2                           GO_CARDIAC_CHAMBER_DEVELOPMENT POLR2J
3 GO_DNA_DEPENDENT_DNA_REPLICATION_MAINTENANCE_OF_FIDELITY  CTDP1
4                                      GO_CIRCADIAN_RHYTHM   RDBP
5              GO_PHOSPHATIDYLSERINE_ACYL_CHAIN_REMODELING COBRA1
6                               GO_SPINAL_CORD_DEVELOPMENT   RSF1

Edited_response
You can use the code given below to achieve your desired output:
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
len_vec=c()           # Now create a vector for containing the length of genes at each position
len_vec[1] = 3
for(i in 1:length(data$genesets)){len_vec[i] <- c(length(data$genesets[[i]]))}
pathway_vec <- unlist(Vectorize(rep.int)(data$geneset.names, len_vec),use.names = FALSE) # Now create a vector for all the pathways in the data
desired_df <- as.data.frame(cbind(pathway_vec,unlist(data$genesets,use.names = FALSE))) # This gives your desired dataframe
head(desired_df)
                                    Pathway      Gene
1 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2J
3 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION  CTDP1
4 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RDBP
5 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION COBRA1
6 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RSF1

sorry, this is wrong. Please check the output I want in the above question. In your output Genes POLR2J, CTDP1, RDBP, COBRA1, RSF1 shows different pathways.


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:
Accessing UCSC genome via ssh results in a validation error
I am working on macOS High Sierra.
I am following the steps described here to lift an annotation over from one version of a genome to another. I am now just using the example genomes provided in the tutorial.
At this step:
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

I get the following error in the do.log file:
Warning: Permanently added 'localhost' (ECDSA) to the list of known hosts.
username@localhost: Permission denied (publickey,password,keyboard-interactive).
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Thu Nov  8 14:27:32 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
align: localhost does not have /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.ooc -- if that is not the correct location, please run again with -ooc.

I think the main problem is here:
username@localhost: Permission denied (publickey,password,keyboard-interactive).

This post describes this issue. However:

I do not have this file: ~/.ssh/authorized_keys. My files in ~/.ssh are: github_rsa github_rsa.pub known_hosts. One of the lines of known_hosts is localhost. Does this presumably allow me to ssh to my own laptop? I still don't understand why I need to ssh in this script if it's run locally.
PasswordAuthentication yes is already in /etc/ssh/sshd_config
If I run the command without the time command, I still get the same error. I am not prompted to give a password.

Could anyone help me overcome this problem, please?

You need to setup ssh keys.  All the commands in the script are going to use ssh to run the commands.  Note the tips here of how to setup ssh keys:
http://genomewiki.ucsc.edu/index.php/Parasol_job_control_system#SSH_keysssh key setup


text:
How to get PSSM of protein by protein sequence with python?
I didn't find way that convert protein sequences in  FASTA to PSSM by python,
and how to write program in python to find PSSM for our data set ?
Maybe better to question how to find PSSM from NCBI alignments ?
Thank you

Welcome to the site. Where did you search? Did you look up to biopython? What is your data set? A collection of fasta sequences or alignments?

Thank you for your answer , yes i check biopython , but there is no method or class that get protein sequences in fasta and return pssm  . for example from NCBI site , for each sequences.

Ok, so you want to retrieve from the NCBI site, download the protein sequences and then create a pssm? That is different from what you posted. Please [edit] the question to include this details. Would you be fine with a solution that worked with a local blast program?

I think better to use psi-blast from NCBI to find my alignments for each sequence  i don't know easy way to  get aligments and pass them to function to get PSSM (i dont know it is necessary to find pssm)

Sorry I don't understand your last comment. What are you trying to do/get? To get a PSSM you need an alignment of sequences. So you need several sequences and then align them (that's what psi-blast does). But you seem to ask about retrieving the PSSM matrix from a psi-blast search. Have you tried something else like hmmer?

Sorry , i am too beginner . Can you give me  a solution to get PSSM for each of my protein sequence with use psi blast NCBI ? I can calculate PSSM localy on my computer but how to get aligmnets?

Sorry it is not easy to answer your question. The alignments must be done downloading all the sequences of interest and aligning them via t-coffee or similar tools. Then you can calculate the PSSM matrix (not the other way round)

@Llopis Thanks, do you know how to automatically download alignments for each protein sequence ?

Download from where? BLAST perhaps? In what format, and in what language ?

From NCBI , PSI-BLAST, in xml fomat, in python

This would require a major re[edit] of the question. Could you please update the question to keep the question and answer related?

We solve our problem and create gitHub repository (bio-protein).


text:
different results coming from biomart online and biomaRt R library
I am playing with biomart in order to get the list of genes associated to a specific gene ontology (GO) term.
So, if I use the webservice as in this query:
I get at least 50 genes. If I download the file I can gen an exact count of them:
$ cat mart_export.txt | wc -l
347

Now, I can do the same operation by using the R library biomaRt. Specifically:
> ensembl <- useEnsembl(biomart="ENSEMBL_MART_ENSEMBL", dataset='hsapiens_gene_ensembl')
> gene <- getBM(attributes = c('external_gene_name'), filters = 'go', values = 'GO:0030098', mart = ensembl)
> unique(gene)
  external_gene_name
1               SPI1
2              PRKDC
3                TOX
4               LY6D
5               CBFB
6              IKZF1
7               RELB
8             FLT3LG

By doing so, these are all the genes I retrieve. The two results are totally different and I don't know why it. All these eight genes are found in the biomart file (I grepped it, to be sure).
I believe this is not an issue related to databases versions since I checked:
>listMarts()
               biomart               version
1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 94
2   ENSEMBL_MART_MOUSE      Mouse strains 94
3     ENSEMBL_MART_SNP  Ensembl Variation 94
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 94

and:
>listDatasets(ensembl)
...
38 Human genes (GRCh38.p12)
...

Both are the same version used by biomart online.
What is going on? Is it possible that biomart online is considering GO relations that are not specified in the R version?

Funny: if I count the result of your query I get that there would be 361 genes. While we have 346 genes (you need to discount the header, and also as a tip wc -l mart_export.txt works as well). The difference between 361 and 346 might be duplicated gene, which suggests that they are considering GO child process. But I think that Emily might have a better answer for this. Have you searched online for similar questions?

yes, you're rigth about the header ;-) I forgot about it. also yes I looked for other similar questions and I found the one that suggested me to check the version used by R biomaRt. Another similar question https://www.biostars.org/p/14948/ does not answer. How come that you get 361 lines? I'll try that again, sending another equal query. Lets see

Ok, I tried again and I still got 346 genes! This is weird!

The 361 is using the count feature of the query page. Instead of results I pressed count at the top left header. If you deselect the unique rows you'll get more genes

yes, I tried it but actually I don't know. Hopefully unique genes will be enough!

Sorry, this was my mistake in the last question. To search down the ontology, rather than just for the specific association with a term, biomaRt needs a different filter: 'go_parent_term'.
Try:
gene <- getBM(attributes = c('external_gene_name'), filters = 'go_parent_term', values = 'GO:0030098', mart = ensembl)

ooook, that works perfectly!
A small point, though. Is this using only "is_a" relations or is it using any kind of GO relations?

These are the mapping we use between parent and child: biological_process = is_a,part_of,regulates,positively_regulates,negatively_regulates;
cellular_component = is_a,part_of;
molecular_function = is_a,part_of


text:
How to link GDC ids to CCLE cell line names?
Hello bioinformaticians,
I've recently downloaded few RNAseq data files from Genomic Data Commons Data (GDC) portal. These files belong to Broad institutes CCLE project.
Now the problem is that GDC portal do not provide CCLE cell line names:

So my question is how could I link the GDC ids (prefixes of File Names) to CCLE cell line names?
I would imagine because it is an important step in the analysis, there should be some way to do that.

Where have you searched? Is there any information page about how to identify the cell lines or the files?

It's an issue with GDC portal. So I looked around on their website. There is not even an email address to query about this issue.

The filenames of those bams contain the cell line names. They are constructed:
SAMPLEID.CELL_LINE_NAME.REPLICATE.bam
Thus the first 5 cell lines in your screen capture are SK-N-F1, TE-9, CJM, HCI-H1915 and CAL-33.
All the information that is need is also present in the experimental metadata xml file attached to each sample record. This is accessible by downloading it from the record on the website, or by querying the GDC metadata API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api).

Thanks. There is indeed some information in the filenames. However, the CCLE cell line names  are generally in the format {CELL_LINE_NAME}_{TISSUE_NAME}. For example, CJM_SKIN and CAL33_UPPER_AERODIGESTIVE_TRACT.
I'm not sure where can I get the tissue names now.

If you look at the "project" column in the above table, it gives the cancer type code from which the cell lines were extracted.  Now how to map those abbreviations is less obious, but if you click a primary site filter in the left hand column, then you can work it out.

This information is also present in the experimental metadata xml file attached to each bam file. This probably means its also accessible from the GDC metadata API.

Ok. I see. There are XML files along with the BAM files.

GDC API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api) should be the best way to dig into those. Thanks.

Could you add the xml file related comment in the answer so I can accept it as the correct answer?

I would recommend using UCSC Xena instead of GDC for CCLE data. They don't have everything. For example, BAM files are not available. However, the data that is available is organized in standard text files with consistent identifiers.


text:
How to run same command on multiple files?
Hi im doing Variant calling on fastq files for which i have 4000 fq files and the variant calling are done in different 9 steps. Each step generate different files that are output to next steps.
So can anybody help me with this to do it in a loop it pick up names from previous steps for all files. Step 1 command is here:
bwa mem -M -R '@RG\tID:ERR040140\tLB:2517707\tPL:ILLUMINA\tPM:HiSeq2000\tSM:ERR040140' ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta ERR040140_* > ERR040140.sam

So i made an excel file to break the command in columns so that i can change the file names keeping remaining command same just like this:

So changes are occurring 5th, 7th and 9th column. which are:
tID:ERR040140
tLB:2517707
tSM:ERR040140
ERR040140_*
ERR040140.sam

Can anyone help me with loop or something?
Thank you!

Please [edit] your question and replace the image of the excel file with a copy/paste of the same file exported to tsv (save as 'tab separated text file'). That way we can actually copy it and use it to test our examples.

Please remember to confirm an answer once you've received one.

I have already given the excel file command above. And im unable to paste it here the page is not letting me paste it that is why i quoted 1 example and than added a picture with all other repeating commands with different file names.

Another solution for repeating one command with different parameters is to use parallel.
Create a tab-delimited file samples.txt, which contains the sample names and the information for your read groups, e.g.:
Sample1 Lib1
Sample2 Lib2
Sample3 Lib3

Then you can use parallel like this:
parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:{1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta {1}_* > {1}.sam' :::: samples.txt

For each line in samples.txt a new process is started, replacing {1} and {2} with the values in column 1 and 2 of this line.
fin swimmer

wow that's look much easier i will try and update you.

swimmer thank you for this simplest code. I want to know whether this parallel produce log? because i cant see anything on the screen as the bwa command produce screen messages while running. Secondly, if i don't want to run these commands parallel than what should i used instead of parallel as while running it parallel my computer got stuck.

By default parallel prints to the screen after a job has finished. Doing this makes sure that you can e.g. concatenate outputs of parallel jobs without mixing them. You can use the option --ungroup to get the output immediately. But use this with care! Without any parameter the number of parallel jobs is the same as the number of CPU cores. You can set it to an explicit number using -j, e.g. -j 1

Ok. Thank you. It is running for this command i will modified it accordingly for other steps and if i got stuck i will update here.

Something like this might do it:
#!/bin/bash

# make an array to hold your list of samples/base filenames
samples=( S1 S2 S3 S4 S5 )

# make a loop to go through the sample IDs
for sample in "${samples[@]}"
do
  :
  # do something with each sample/file, e.g. make a directory for it
  mkdir "${sample}"_home

  # do something else with the sample, e.g. unzip the fastq.gz file
  gunzip "${sample}".fq.gz > "${sample}"_home/"${sample}".fq
done

One disadvantage with this is that it will process the files one after the other.  If you have a cluster, you might prefer to run them all at the same time to make the whole thing quicker.  If you have a SLURM cluster, you can use job array to do something similar.
If you have your list of samples in a text file, you could try replacing the samples array command in the script above with something like this:
samples=( `cat samples_list.txt` )

This will only work if the samples are listed on a single line, space-separated. You may need to process your samples file further to get them in this format, or use the bash readarray command.

I have a list of sample ids store in a txt so is it possible that it make a list using that file?

@SafinaA.R please remember to actually show your data. Edit your question and show us this text file. That way, we will know the format and can use it to give you example answers.

This solution requires saving the sample names directly into the bash script. This sounds a little difficult since OP said she has a file with 4000 sample names.

hi @terdon
The data is same as shown by fin swimmer in the above example. It is a text file tab separated having two col. Sample Ids and Library Ids

I haven't tried this, but I added a comment to my script above about how to get the sample IDs out of a text file.

if you have a repeating workflow, I strongly recommend to have a look at workflow management systems like snakemake. I've also wrote a little tutorial on biostars about this topic, which might be useful for you.
fin swimmer

That's a bit complicated can there be simple for loop trick?

You can create a script to pass arguments to as follows:
#! /bin/bash

bwa mem -M -R "@RG\tID:${1}\tLB:${2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:${1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta ${1}_* > ${1}.sam

Now if you run bash script.sh sample lib then sample and lib will be  the script as arguments ${1} and ${2} respectively. Then all you need is a loop to pass arguments to the script:
for sample in sample1 sample2 sample3
   do
   for lib in lib1 lib2 lib3
     do
     bash script.sh $sample $lib
     done
   done

You can use nohup to run each sample in the background and proceed with next one. This will in essence, submit each command as a separate job to run on a different core. These will print to separate log files sample.out instead of to the terminal. It will also continue to run, even if you disconnect from SSH or close the terminal.
for sample in sample1 sample2 sample3
   do
   for lib in lib1 lib2 lib3
     do
     nohup bash script.sh $sample $lib > ${sample}.out &
     done
   done

You should only do this for single-threaded processes. For multi-threaded processes or memory-intensive jobs, this will demand a lot from the machine. Even if this isn’t a shared server, it will perform poorly if you run out of memory and have to “swap” memory by writing to disk. So it’s best to submit them one-by-one (above) or use a cluster with a job scheduler.


text:
Bismark report before and after Cytosine Methylation Extraction
Looking at the sample report from Bismark webpage, there are two sections: Cytosine Methylation and Cytosine Methylation after Extraction. However even after reading User Guide, I still cannot understand the difference between these two sections.
As far as I understand this, Cytosine Methylation section shows summary statistics about Cytosine Methylation in sequenced data and Cytosine Methylation after Exctraction shows summary statistics after a Methylation calling, but is not the Methylation calling already part of the first section?
My question therefore is: What is the difference between these two sections of Bismark report?

Did you try contacting the authors of Bismark? It doesn't seem to be highly used, so I doubt we have here someone who could answer you. (If you find the answer, then post it here for others please)

I was also looking for an answer to the same question and found this issue. have a look at this post https://github.com/FelixKrueger/Bismark/issues/245


text:
How to plot a 2D RNA structure in python using dot-bracket input?
I need some help with forgi library for visualisation of RNA secondary structure.
I have an RNA sequence and dot-bracket notation of the secondary structure and I would like to plot 2D graph like that you can find here.
The input file is a .fx file, or a fasta file where the sequence is single-line and the third line is the dot-bracket notation of the RNA structure.
>1y26
CGCUUCAUAUAAUCCUAAUGAUAUGGUUUGGGAGUUUCUACCAAGAGCCUUAAACUCUUGAUUAUGAAGUG
((((((((((..(((((((.......)))))))......).((((((.......))))))..)))))))))

If there any other possibilities to plot it in real-time (for instance, in jupyter notebook avoiding any crossing external scripts), I will be very thankful to learn them!

Hi there, do I understand right that you got an RNA sequence and you with to predict the secondary structure and plot it within a jupyter notebook? Could you maybe [edit] your question to clarify what is exactly your input and desired output?

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

@KamilSJaron, hi!
I have an RNA sequence with it dot-bracket notation, for instance:
seq = 'AAAACCGGGCCUUUUACCCCAAAUUGGAA'
dot_bracket = '((((..(((..)))..))))...((..))'
And for this seq and dot_bracket I want to plot 2D graph where edges will be biult according to dot_bracket and top will show the necessary nucleotide

@kitsune_breeze Ah, I got it. Just note that you can [edit] your question to improve and I would strongly encourage you to do that if you find new clues or to add more details. For instance, what was the problem with forgi library? Could you share what you tried and what was the error you got?

@KamilSJaron, thanks. I cannot really understand how to rewrite, for instance, this command 'python examples/rnaConvert.py examples/input/1y26_ss.dotbracket -T neato | neato -Tpng -o 1y26_neato.png' into something that I can use in jupyter notebook.
I mean what functions I should call from forgi library that allows me to draw appropriate graph

Awesome, that's something we can work at. Just [edit] your question and add this information. Something like, "I know how to create the plot in terminal by running .... However, I don't know how to call the command from jupyter notebook". I am trying to convince you here, that you should not write more information in comments but directly to question. ;-)

I wrote a function to do that a while ago, maybe it helps. It only works with dot-bracket notation, and it turns it into ordered pairs which then you can easily plot.
import pandas as pd
import seaborn as sns

def dot_bracket_to_pairs(ss_string):
    '''
    Takes dot and bracket string, returns dataframe
    with paired bases.
    If any invalid characters are in the structure, it
    will interpret them as dots, as it only reads parentheses.
    '''
    index_list = []
    pairs = {}

    for index, char in enumerate(ss_string):
        if char == '(':
            index_list.append(index)
        if char == ')':
            try:
                # pair to last item in the list in dictionary
                pairs[index_list.pop()] = index
            except IndexError:
                print(f'Invalid structure, found extra \')\' in position {index}')

    if len(index_list) != 0:
        for item in index_list:
            print(f'Invalid structure, found extra in \'(\' in position {item}')

    df_pairs = pd.DataFrame(pairs.items(), columns=['Y', 'X'])

    return df_pairs

db_string = '((((((((((..(((((((.......)))))))......).((((((.......))))))..)))))))))'
df = dot_bracket_to_pairs(db_string)

sns.scatterplot(data=df, x="X", y="Y")

A cool, yet simple, solution, thank you. The plot is y-axis 5' and x-axis 3' ... switching would be clearer. A lot of coding is just the concept and this is cool.

One way simply is to use VARNA, here is a link ... https://varna.lri.fr. This is well known for its power to view secondary structure (e.g. its in Jalview). It is generally considered the best viewer app for secondary structure.
However, the solution by @nquinoneso whilst simple, is preferable. The output is here:

The y-axis is the 3' and the x-axis is the 5', the dots represent hairpins. Whilst its simple, the power of this approach is the automation by direct access to the (pandas) dataframe.
I do recommend upvoting the post.
Complex secondary structure would not work because it will confuse the 5' stems with the immediate 3' stems which don't hybridise and  a more advanced solution using nucleotide code would be needed. Thats the bit for further development, using the index both for the dot-bracket and the nucleotide sequence to confirm the hairpin. VARNA does perform this calculation (because it identifies canonical mutations)... but you can't access the underlying array/dataframe which is certainly a limitation.

Sorry, I got the x,y the wrong way round. The y-axis is the 5' end and the x-axis is the 3' end ... making the plot difficult to read. Obviously easy to correct.


text:
Mapping Reads to Known Gene Paralogs with Long Read Technology
I have some sequencing data from a captured region that is a known paralog edited. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:
minimap2 -ax map-pb GRCh37.fa fastqfile.fastq > samfile.sam

I captured a 10kb region and what I am seeing in IGV is puzzling. The region directly before and directly after the paralogues are covered at about 1000X, but the genes themselves are not covered at all:

The two genes in the middle are the paralogues. I left them off to protect data. This data is also only the sorted data, the duplicate marked data has lower coverage (600X).
Question 1
I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?
Question 2
Is this process normally accomplished by minimap2, or are other aligners more effective for paralogous mapping?
Question 3
The interesting thing is, while my average read length was 6kb, I am seeing the reads here are about 4kb. I have another region with very high coverage where the reads spanned both paralogues and mapping was successful.
Question: Why am I only seeing reads that are about 4kb long when the average read length for this dataset is 6kb?
Any help on this would be great!
Thanks and I hope I have been clear on my question,
Dennis

Given that there are no spanning reads, I think that your igv screenshot is more likely to be caused by non-specific capture than by alignment. If it were me, I would double check capture probes.

Question 1
This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.
samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam

This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the bam file.  Using picard tools is the fastest way I have seen to filter a bam file given a list of reads.
#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$filtint" ${BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq > ${TEMP}; \
java -jar ${PICARD_JAR} FilterSamReads \
  I=${BAMFILE} \
  O=${OUTFILE} \
  READ_LIST_FILE=${TEMP} \
  FILTER=includeReadList
rm ${TEMP}

This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.
Question 2
Filtering the bam file should be pretty standard, as addressed above in my response to Q1. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.
The fact that there are paralogs does not matter to the aligners. The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. I do not know of any long-read paralog-aware aligners.
As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.
Question 3
Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.
One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.

I will definitely give this a try! Thanks for the input. I will upvote this for sure but am looking for a more complete answer for handling these paralogues. I am giving blasr a try and trying some titrations of input variables. It has some advanced parameters for largely identical sequences (which is what these are). Thanks!

Thanks. If you are able to edit your question to be more clear about how you would like to handle the paralogs we can try again. Right now from your question this appears to be the only sentence that says what output you would like: "I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?"

Sorry for the ambiguity, I have updated the question to address the idea of different mappers or approaches. Thanks

Great, complete answer. Thanks for the edits and update!

I have found an alternative solution that has been excellent for this paralog mapping issue, so I thought I'd post what I had worked on:
Question 1+2
I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?
Using blasr, pacbio specified some options in their github docs for long, nearly identical reads:
blasr input.bam /path/to/reference.fa --hitpolicy all --bam --out alignments.bam --minMatch 8 --maxMatch 15 --nproc 4

The link to the docs as to why this is useful is here, although I did tweak the parameters to smaller seed sizes to create more anchors to the reference, which produced more secondary alignments. This is effective because this will output up to 10 secondary alignments per read, which I saw quite frequently for my region of interest.
After this, sorting was accomplished using Picard SortSam, with the following flags:
java -jar picard.jar SortSam INPUT=${SAMPLE} OUTPUT=./sorted/${SAMPLE}\
SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT\
CREATE_INDEX=true MAX_RECORDS_IN_RAM=100000

The blasr files have much more information than normal sams produced by aligners like minimap2 and bwa mem, so the MAX_RECORDS_IN_RAM=100000 as opposed to 500000 helps to ensure less frequent job failures (this was done on a cluster).
Next, I parsed out the region of interest:
samtools view -h input.bam "chrX:1000-2000" > extracted.bam # this is fake

Parse out mapq=0 reads:
awk '$5==0' extracted.bam > mapq0.txt

Get read names to grep for later:
awk '{print $1}' mapq0.txt > read_names_mapq0.txt

Search for all reads with readname from the original bam to get mapping positions:
samtools view infile.bam | grep -f read_names_mapq0.txt > all_reads_mapq0_mapq_high.txt

Keep only columns: 1) read_name, 2) Sam Flag, 3) Chr, 4) Pos, 5) MapQ:
awk '{print $1,$2,$3,$4,$5}' OFS='\t' all_reads_mapq0_mapq_high.txt > mapq_read_info.txt

Then I wrote a simple python script to sort by read name, count the number of times a read shows up, and see all the positions it maps to. An example of the output:
        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       272     19      3395674 0
read1/10158955/ccs       0       16      217682  254
read1/10158955/ccs       272     12      64990527        0
read1/10158955/ccs       256     16      222913  0
read1/10158955/ccs       272     19      21328060        0
read1/10158955/ccs       256     19      3224052 0
read1/10158955/ccs       272     17      34878621        0
read1/10158955/ccs       272     16      84204312        0
read1/10158955/ccs       256     1       244788680       0
read1/10158955/ccs       272     2       222951342       0
read2/10223868/ccs       16      16      217057  254
read2/10223868/ccs       256     20      3078724 0
read2/10223868/ccs       272     6       133521134       0
read2/10223868/ccs       256     7       152515722       0
read2/10223868/ccs       256     21      45258117        0
read2/10223868/ccs       256     1       249160685       0
read2/10223868/ccs       272     16      223631  0
read3/10551872/ccs       272     GL000220.1      144982  0
read3/10551872/ccs       256     18      72831954        0
...

read_counts output:
        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       10      10      10      10
read2/10223868/ccs       7       7       7       7
read3/10551872/ccs       9       9       9       9

The numbers in each column correspond to the number of supplementary reads seen.
A brief explanation of sam flags which I feel is relevant for unpaired reads:
Flag=0: read was mapped to + strand
Flag=16: read was mapped to - strand
Flag=256: read was mapped to + as secondary alignment
Flag=272: read was mapped to - as secondary alignment

Minimap2 gives you up to 5 secondary alignments by default. IGV just doesn't show them.

Yeah I know, but I tried many varying -H and -k parameters and never saw any secondary reads to my region. I know they should be there because when I chunk my reads and blast them I see 100% identity to this region as well as a few others but minimap2 is not outputting any secondary alignments there.

And I don’t just mean in the viewer, I mean in the sam file. Unless there is more I can do?

100% identity doesn't count if the alignment only covers a small fraction of the query. Those fragmented alignments are probably hits to repeats. Minimap2 can be tuned to output such hits with -p.1 or -P, but most of extra alignments are not useful.


text:
How can I identify a recessive and dominant gene?
For example, I have two allelic genes. How can I identify a recessive and dominant gene? Are there any databases with this information? Or the answer to my question is to study the concept of "haplosufficiency"?

I'd probably look up the genes in OMIM, which provides detailed description of the functions and phenotypes associated with genes, including inheritance mode.
For example, this is what you would get for TP53:


text:
mrbayes, Sumt ; Error when setting parameter "Filename" (1)
question 1
I wanted to use sumt command after completion of mcmc run. Whenever I use, sumt Filename=infile.nex, it gives me the error:
A matrix must be specified before sumt can be used
Error when setting parameter "Filename" (1)

But for the same analysis, I can use sump without any problem. What should I do now?
question 2
Another issue, I used burn in as 25% where as in Tracer, it is always showing 10%. No matter how much burn in I use, Tracer always shows 10% as burn in.
system stats
Operating system : Centos7
Version of MrBayes: 3.2.6

Thanks in advance,
Best Regards,
Zillur

Hello zillur rahman,
you've crossposted this question also on biostars. You should mention this, so that everyone can see, if there is already an answers to your problem.
Thanks.
fin swimmer

Having thought this through, I think you've made a mistake. The file it wants is the outfile, you have given it the infile only. My original answer is below.
====================
From memory (my personal memory, not RAM, SSD or Stata), in the run command, you need to specify,
nruns=2

sumt will then combine both runs. If the file name is mystuff.run and nruns=2, there will be an output of mystuff.run1.t and mystuff.run2.t.
If nruns isn't set then sumt will not work, or if the MCMCMC output ends up with different names again sumt will not work.
With the burnin, you simply specify this in MrBayes. It is easy to write a script that will perform a 25% burnin. Tracer is nicely written, presumably C, its fast, but it has been around for a long time.

This is a bit late. However, I am hoping someone else will find this useful since I had the same problem.
The data matrix is your infile.nex file (the nex file you supplied for the analysis)

Bracket out the mcmc command in your infile.nex

Prevents the file from being run again.
May not be present.

In Mr. Bayes: Execute infile.nex

Loads you infile as the data matrix

In Mr. Bayes: sumt Filename=infile.nex nruns=2

If you have nruns set in you Mr. Bayes block you can skip this.

This question was posted and answered as a Github Issue in the MrBayes repository; I am copying the answer here for better visibility.

Dear zillurbmb51,
Thank you for your question. You need to load the data matrix in
memory before issuing the sumt command, or in other words: "A matrix
must be specified before sumt can be used".
Example 1 (where you haven't yet run a MCMC):
Start mb, then type:
exe infile.nex
mcmc
sumt

Example 2 (where you already did a MCMC, and have .t files in your
current working directory, and they are named infile.nex.run1.t,
infile.nex.run2.t, etc):
Start mb, then type:
exe infile.nex
sumt filename=infile.nex

Cheers
Johan


text:
Why is a PacBio read length larger than the aligned reference region?
I recently had some Iso-Seq sequencing done on my organism catfish on the new Sequel platform and got weird alignments for a size selected 4 Kilobase and up fraction after running the isoseq3 pipeline.  When I investigated the exon intron junctions using data from (Pac Bio Iso-Seq, Illumina RNA-seq) against our private reference genome, I get PacBio Iso-Seq reads where the sequence length is greater than the seq length on the reference for my bam alignments:
Seq length: 6751
Seq length on ref: 5267
When blasting the sequence in NCBI it shows 5% similarity to the publicly-available genome.  I'll attach a JBrowse image with the same RNA sequenced with Illumina paired RNA-seq.  I also included a image and circled the pacbio contig and Illumina is above.

Question 1:
Is this a problem with JBrowse?
Question 2:
Does this mean that some DNA was sequenced with the catfish cDNA?
Question 3:
Why is this aligned sequence's length longer than the length of the reference genome region to which it is aligned?
Thanks

Q1
Probably not - this is probably just displaying the contents of the input (I presume a bam file) exactly as they are encoded.
Q2
The splicing looks a little weird, but the exon-mapping pattern in the 5' part of the gene makes it seem like RNA rather than DNA. I think that the large intron is what makes you think that it is DNA, correct?

Explanations for weird splicing

The weird read is a rare splice variant but a true mRNA.
The weird read is not completely matured and has an intron that has not yet been spliced. The IsoSeq library prep protocol picked this up somehow and it ended up being sequenced.

Q3

Explanations for why Seq length: 6751 and Seq length on ref: 5267

There is probably a large amount of clipping on the 5' or 3' end of the read. This could be polyA or some other transcribed-but-not-yet-completely-matured seq. Or, the extra sequence could be sequencing error, PacBio adapters, or a chimeric molecule.

things to try
How to check what is going on:
I'm guessing that you had to upload a bam file to JBrowse to make the IsoSeq track. Use the sequence name that you got from Jbrowse and run the following command to extract the alignment info from the bam file:
samtools view /path/to/your/isoseq/file.bam | grep 'weird_read_id'

Look at the cigar string to see what it says about 5' or 3' clipping and corroborate with the sequence data.

you are right the first 3702 bases are soft clipped.

I'm curious to know what the clipped sequence is once you figure it out!

I am still a mystified why it is clipping.  I am starting to notice that many of the reads that display this property look chimeric where the start site is seen twice in the middle of the sequence. I also note that this was the 4kb and up fraction.  So it might be assembly error you mentioned earlier I am feeling.

Thanks for following up. Please feel free to post another question about the assembly error if you'd like!


text:
Why are TPMs per 10k or 100k in many scRNA-seq studies?
I noticed that many scRNA-seq papers normalize TPMs to 10k or 100k as opposed to 1M (as the abbreviation defines them). It doesn't really matter since you are just moving the decimal point, so why mess with an established convention?
Additionally, the values might be log-transformed. In that case, log(x+1) would affect the data less if the total is 1M.

Could you provide some links, so that we can check if there is some explanation in the methods section?

@GWW Isn't "massive" subjective? With standard TPM, the highest possible value is 1M. It doesn't seem that big. Arguably, by switching to 10k, your small values become too small (extra 0s at the beginning).

When studying demographics, you often need to account for the population size of a city or region (GDP per capita, is a good example).
Normalizing by $10^6$ people is a nice round number, but lots of cities and regions don't have that many people.
So if we want to keep numbers in the range [1, 100], 1M is too big a denominator.
That's why a lot of studies use $10^5$ people as a scaling factor.
The same goes for RNA-seq and other sequencing measurements.
Your readout will typically contain $10^7 - 10^8$ reads for a bulk sequencing run.
If your sample has thousands of cells, this means very few cells will have > $10^6$ reads assigned to it.
But the TPM unit was designed for bulk measurements and that's what most people are used to.
So we can keep the same idea of scaling by number of reads, but just use a different convenient scaling factor.
This way we get both scaling by sequencing depth and nice numbers to work with.
As for your $\log(1+x)$ comment, yes, if $x$ is in TP10K, the low read counts will be biased more than if $x$ was in units of TPM.
But consistently lowly-expressed genes often aren't of interest.
So this won't be too much of a problem.
The bigger issue is the zero count over-inflation that is common in scRNA-seq measurements.

Collecting all the question comments as a community answer (feel free to edit this to make it more readable):
There tend to be far fewer counts in single cell experiments. By reducing the scaling factor you are avoiding having massive expression values for your samples.
It isn't uncommon to scale single cell experiments by the median cell library size. That way small values typically won't be too small and large values won't be too large. People also choose 10,000 because that's roughly around the library size for a single cell.
I haven't seen an explanation yet. Here is an example paper from last week (reputable journal and lab): doi.org/10.1016/j.celrep.2018.11.056 . "... We also excluded outlier cells from the MPPs group which were likely caused by cell contamination. To correct for batch issues, the two experimental batches were centered to have the same mean log(TP100K+1) (hereafter referred to as log(TPM), where TPM stands for ‘transcripts per million’) per gene. PCA was done in R with scaled..."


text:
Cicros plot issue using two data track
Im trying to plot log fold change for a condition as well as no of patient for a particular gene mutated information but whenever im trying to do that i see only only the fold changes labeling not the patient but my idea is to see the fold change as well as mutation information for that gene im posting the .conf file along with the data files im using any help would be highly appreciated .
Or there is a flaw in my idea to show both of them which im not sure but i guess whatever i seen so far examples it can work!
My conf file code and my data files data
use_rules = yes
karyotype = data/karyotype/karyotype.human.txt
# spacing = 10000000
#
# you can write
#
# spacing = 10u
#
# See ticks.conf for examples.

chromosomes_units = 1000000
#chromosomes_display_default = no
#chromosomes                 = /hs[1234]$/
#chromosomes=hs5;hs7; hs17; hs11; hs9

<<include ideogram.conf>>

<<include four_ticks_4chr.conf>>

<image>
<<include etc/image.conf>>
</image>
<<include etc/colors.conf>>

<<include etc/colors_fonts_patterns.conf>>

<<include etc/housekeeping.conf>>

<plots>
<plot>
#this is for two direction expression data visualisation
type = histogram
#file=data2/HSC_VS_CMP.txt
#file = /home/punit/circos-0.69-5/data_files/LEUKEMIA_value_tab.txt
file = /home/punit/circos-0.69-5/data_files/CIRCOS_VAL1.txt

r0   = 0.60r
r1   = 0.35r
thickness = 10
min=-1.5
max=1.5
extend_bin = no

<rules>
<rule>
condition    = var(value) > 1.5
color        = dred
fill_color   = dred
</rule>

<rule>
condition    = var(value) < -1.5
color        = dgreen
fill_color   = dgreen
</rule>
</rules>
</plot>
#########################################################

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

<plot>
type    = text
color = black
#file = /home/punit/circos-0.69-5/LEUKEMIA_gene_tab.txt
file = /home/punit/circos-0.69-5/CIRCOS_LOC1.txt

r0   = 0.75r
r1   = 0.95r
show_links     = yes
link_dims      = 0p,0p,50p,0p,10p
# link_dims      = 4p,4p,8p,4p,4p
link_thickness = 2p
link_color     = grey
label_snuggle             = yes
max_snuggle_distance      = 2r
label_size   = 25p
label_font   = condensed
padding  = 0p
rpadding = 0p
</plot>

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

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

    <plot>
    #this is for two direction expression data visualisation
    type = histogram
    color = blue
    #file=data2/HSC_VS_CMP.txt
    #file = /home/punit/circos-0.69-5/data_files/LEUKEMIA_value_tab.txt
    file = /home/punit/circos-0.69-5/M_DATA.txt

    r0 = 2.0r
    r1 = 1.8r
    orientation = out

    #max       = 20
    #min       = -20
    stroke_type = outline
    thickness   = 18
    extend_bin  = yes
    fill_color = blue

    <rules>
    <rule>
    condition    = var(value) > 10
    color        = dblue
    fill_color   = dblue
    </rule>

    <rule>
    condition    = var(value) < 8
    color        = dyellow
    fill_color   = dyellow
    </rule>
    </rules>

    </plot>

    </plots>


text:
scRNA-seq, 10x cellranger pipelines
I'm starting to do sc-rnaseq using 10x cellranger pipelines, and i add TdTomato sequence to mouse reference genome and add an entry in the gtf. My code is But when I makref, it remindered that:
Writing genes GTF file into reference folder...
Invalid number of columns in GTF line
Can anyone tell me the reason?

Do you have a fasta and gtf file for your custom gene? It sounds like you can just add them below your genome.fa and genes.gtf files respectively with cat command.

I have add  fasta and gtf file to  genome.fa and genes.gtf files, but cellranger did"t work ,and I update my problem and related pictures,can you heip me find the reason?

The error messages you shared within your question and within what you have pasted is different. The first error is about line 1794438, which is 20 . exon 1 1431 . - . gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato". As you can see, the first field of that line is not the same as your gene name that was used in the fasta file, but is "20". I guess it is also a good idea two double check the spacing between different fields in your GTF, which is causing the problem and should be tab-separeted.

yes,thank you.  I should check the spacing between different fields in  GTF.but the first field of 1794438 line is 20,( 20 . exon 1 1431 . - . gene_id "TdTomato"; transcript_id "TdTomato"; )why you say is "20", I did't find.

thank you all for your reply and help, the reason for the bug is when I copy gtf code from windows to linux, the tab changes,and gtf was not tab-separeted.Thank you again!

@sophia,I see you replaced your initial question and made it into a new one, Please leave the original question intact, and start a new question for your other problem.

Please don't edit the question so that it asks something different! That invalidates the existing answers. If you need help with the step, post a new question. And please never post screenshots of text like that! We need to be able to copy your data and code so showing us an image is essentially useless.

As was suggested by @benn, you will need to add the TdTomato sequence to the end of your fasta formatted genome file. If your TdTomato sequence is in fasta format already, cat will do the trick. If not, just go the end of the file, start a new line and add a header line starting with ">" then immediately the name of your gene, no whitespcaes in between, something like ">TdTomato". Then add your sequence in a new line.
The format of the GTF file is explained in the image you posted, the trick is to use the exact same name (in column 1) that is used in the genome file you will create at the previous step (no ">" though). Again, this line, in which the different fields (columns) will be separated by tabs, will be added to the end of your GTF file.

Thank you for your detailed answer,  and I have add  fasta and gtf file to  genome.fa and genes.gtf files, but cellranger did't work ,and I update my problem and related pictures,can you heip me find the reason?

thank you all for your reply and help, the reason for the bug is when I copy gtf code from windows to linux, the tab changes,and gtf was not tab-separeted.Thank you again!


text:
Identifying Indels from Chromatograms
I have around 100 chromatograms (.ab1 files) from Sanger sequencing a genome at loci believed to have an indel.
I'm new to interpreting this kind of data in general, but I've read a bit on the general idea—mostly in guides like this. What I'm not able to figure out from the resources I've checked is how to identify insertions and deletions in a chromatogram.
I'm most concerned with  heterozygous indels, and these seem to be simpler. (And since my samples are (supposed to be) larger indels that are more rare and more likely to be het, I think.) A heterozygous indel, according to my understanding, would like something like consecutive heterozygous SNPs—overlapping peaks, maybe not with the same height, but with the same horizontal placement—until the end of the sample.
In other words, something like the top track in the open window in this image.

(Source: CodonCode.)
(I'm trying to do this manually, and not with packages like this for several reasons, including: I want to see how this is actually done, I want to do it as thoroughly as possible—I will probably check my results with an automated package—and I was asked to do it like this.)
How can evidence of (het) indels like this be distinguished from noise in the chromatogram? My guess is it includes things like

the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and
the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Is this correct?

With my own chromatograms, I have a lot of doubt about distinguishing noise from indel. Here's an example, visualized with 4Peaks. (Full chromatogram here.)

On the left, the chromatogram appears high quality. There are sharp, regularly spaced, single peaks. On the right, there appear to be two sets of peaks lined up with each other. Does this constitute evidence for an indel? Does the homopolymer region—I've noticed in many of these samples, the putative indel seems to start in a homopolymer region—cast doubt on this?

I've also noticed that in some situations like this, there appear to be not two but three peaks.

Is one of these just noise and can be disregarded? Or does this cast doubt on the validity of the data—which is supposed to be coming from the diploid human genome?

If you want to call heterozygous SNPs and INDELs, you may try polyphred. I believe there are a couple of other similar tools, but I have forgotten their names.

How can evidence of (het) indels like this be distinguished from noise in the chromatogram?

Your two points are correct. Basically, at some position, you start to see two streaks of overlapping peaks. Once identified, the two streaks can be aligned to each other if you shift one of them. This is a strong constraint, but it is hard to tell by eye because it is not obvious how to assign peaks to streaks without running an automatic tool. Doing this is even more difficult given noisy traces like your last figure. I think you'd better run polyphred or similar.

the putative indel seems to start in a homopolymer region—cast doubt on this?

The INDEL mutation rate in long homopolymer runs is orders of magnitude higher than elsewhere in the genome. This has two implications: 1) you expect to see more homopolymer indels; 2) PCR polymerase is more likely to introduce errors in homopolymers––a higher fraction of homopolymer indels tend to be wrong.

Thanks! I'll look into polyphred and other tools. To be clear, what's your take on the indel in the homopolymer region—that the chromatogram is inconclusive itself? What can I do to verify it?

@Randoms Overlapping peaks alone are not strong evidence. You need to check if the overlapping peaks can be separated into two sets and the two sets can be aligned with a shift. I can't do this in my head. Even if there is an indel there, the chance of seeing a PCR error is still high. Run polyphred first. PS: finswimmer is right that overlapping peaks should be twice as "low". This makes the second figure less likely to be true.

I think you're right that running a tool is really the way to go. Thanks! I'm going to try Poly Peak Parser: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4525701/

to understand how to interpret sanger sequencing chromatograms I'd like to give a brief introduction to how this data are generated (sorry, if you already know this)
The target region is amplified by PCR. The following sequencing reaction is quite similar to a PCR, despite two important differences:

Only one strand is used as a template, by adding just 1 primer instead of two.
In addition two normal dNTPs, we have modified ddNTPS. Those are fluorescence labeled and the 3'-OH is missing, so there is no chance to add more dNTPS.

The result of this reaction are single stranded DNA fragments of different length. The last base is fluorescence labeled. So we have fragments that looks like this:
A
AT
ATC
ATCG
ATCGA

Of course there a many copies of these fragments.
Now these fragments are separated by capillary electrophoresis. Small fragments are detected first. We can differ the several bases by its fluorescence color.
Let's imagine a DNA strand like this:
ATCGA

What happens if we had on the second strand a SNP?
ATCGA
ATAGA

Now the sequence would detect two signal of equal intensity on base 3 - one for C and one for A. This have two effects in the chromatogram.

Two peaks of equal heights with different colors
The relative peak height is just half the height, compared to those sequences that have just one base at this position.

Let's see what happen in case of a heterozygous deletion. An alignment would look like this
ATCGA
AT-GA

But our capillary electrophoresis doesn't know about alignment. It only knows about fragment length. So this is what it will detect
ATCGA
ATGA

With the beginning of the deletion it will detect to different signals. Just like in case of the SNP, the relativ signal intensity is now just about the half.
With this knowledge let's go back to your initial questions:

How can evidence of (het) indels like this be distinguished from noise
  in the chromatogram? My guess is it includes things like

the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and

They don't start at the beginning, if your indel is not right after the primer, correct. In case of an insertion or deletion you have strands of different size. So if your target region is quite small you may have luck, that you can see the end of the overlap at the end of the sequence, because the shorter strand reaches its end.

the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Yes. And one of the best things to distinguish between noise and true heterozygous variant is a relativ drop off of the signal intensities by approximately 50% compared to homozygous samples.
Concerning your second chromatogram: This is not an indel. As the other already said, this is a mistake made by the polymerase. Invitro polymerases doesn't like homopolymer stretches. They begin to add extra bases or "forget" one. As we do PCR this error accumulates. In my experience up to 9 time the same base is not a big problem, 10 only if you have luck, from 11 upwards your sequence are terrible. If you are interested in the sequence part behind this stretch, you have to sequence from the other direction.
fin swimmer

Thanks for your thoroughness! Interesting that you bring up the relative intensity of heterozygous peaks relative to a het peak at the same spot--I'm not quite sure I'll be able to use this, though, because I've only sequenced one sample, so I don't have multiple sequencing results that cover the same locus.
I'm interested in to know how you're confident that the chromatogram with the homopolymer does not show an indel—surely it's possible that there is an indel in a homopolymer region, right?

You see no significant drop-off of the average peak height and the additional base peaks are to small. That's why it's unlikely that there is a heterozygous indel. But as I wrote, the homopolymer stretch is to long to be sure. I wouldn't give any definitively answer here.

Downvote without comment? Not a good style :(

Sanger sequencing produces a bit of messy signal at the start and a lot at the end of the sequence, but should have a region of clean signal in the middle. If the signal is getting mixed up in the middle, it suggests that there are sequences that are out of phase with each other. That could be a heterozygous INDEL, but it could also indicate that multiple template sequences have been amplified by the same primer.
Sanger sequencing is frequently carried out in both directions (forward and reverse complement), and if similar things are happening in both directions at the same place, it makes it more likely that there is a heterozygous INDEL and not something else funny going on. It's a very good idea to have homozygous templates for all expected variants, if possible. If the secondary peaks in the signal match what is observed in the homozygous variants, that can also be used as evidence. In the absence of control sequences, it is also possible to eyeball the peaks to find self-similarity.
Where triple peaks are seen, that could point to a region where there are two INDEL variants in close proximity, or multiple regions getting amplified at the same time, or even a single INDEL variant with more than two variant states (e.g. a homopolymer with 8 As that is sometimes 7 As, and sometimes 9 As).

Thanks for your response! Can you expand on what you mean by "positive and negative homozygous variant templates"?

Sorry, I was thinking of positive and negative controls; I've changed it to clarify as "homozygous templates for all expected variants"

Honestly I would not try to call heterozygous indels with Sanger sequencing.  It's just not the right technology.
If you were trying to confirm a known indel that didn't flank a homopolymer stretch, and you had sequence from the other side that got messy at exactly the same point, maybe you could safely conclude there was something going on there, and possibly you could eyeball the exact right sequence of mixed peaks, but I don't see that what's what you have in your examples at all.  It just looks like noise.

Thanks for the input. By the other side do you mean the other strand—like forward and reverse Sanger sequencing—or the same strand, just on the other side of the variant?

For me sanger sequencing is still the gold standard for calling indels :) Compared to short-read sequencing is it quite easy to discover indels of any length.

But mixed indels?  Unambiguously no matter what the surrounding sequence?  If there is an indel in that poly-A region, It would be exceedingly difficult to prove.

There are specialized tools to deconvolute (separate) heterozygous InDels in a chromatogram trace using a reference genome. For instance, Indigo is a web application for this, the command-line application is called Tracy.

One important point that @gringer briefly touched on is that there is often off-target amplification when performing PCR on whole-genome DNA.
So, the noise in the chromatograms very well could be from another part of the genome and not related to indels at all.
The only way to ensure that this is not a problem is to clone your PCR products into a vector, transform into something like E. coli, grow up a single colony, cut the product with restriction enzymes, size-select, then sequence the correct-sized insert. Or you could just re-sequence your PCR product with Illumina or PacBio to accurately sequence the complete population of molecules -- indels and all. Check out:
Hebert, Paul DN, et al. "A Sequel to Sanger: amplicon sequencing that scales." BMC genomics 19.1 (2018): 219.
The safe thing to do is to not read too much into the noise of Sanger chromatograms unless it is a simple and very clean SNP.

Before making such a big stunt with cloning, NGS, PacBio etc. one should check the  specificity of the primer, make changes to the PCR protocol, make sure that the clean up steps work, ... If you are working on region that have no homologues, off-target amplification isn't normal.


text:
Truncating branch length values of Phylogenetic tree with biopython
I have been using biopython 1.72 to display my phylogenetic tree files.

Using the function 'Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)' to display branch lengths as well on tree, the tree displays the branch lengths as including all numbers after decimal. While I want the numbers to be truncated to two decimal places.
Following is a python method that truncates upto "the number required" decimal places.
def truncate(n, decimals=0):
    multiplier = 10 ** decimals
    return int(n * multiplier) / multiplier
truncate (12.567, 2)  #Using function
12.56                 #Output
Can I use such method with biopython for branch lengths?
'Phylo.draw' function opens the window which is very narrow and if the tree file is larger, all the branches are overlapped and it gets difficult to read what is written. Does anyone know any better display method?

Following is the example code, I have been using:
from Bio import Phylo
from Bio.Phylo.TreeConstruction import *
from Bio import AlignIO

aln = AlignIO.read(open('example.phy'), 'phylip')

calculator = DistanceCalculator()
dm = calculator.get_distance(aln)

constructor = DistanceTreeConstructor()
njtree = constructor.nj(dm)
starting_tree = njtree
scorer = ParsimonyScorer()
searcher = NNITreeSearcher(scorer)
constructor = ParsimonyTreeConstructor(searcher, starting_tree)
pars_tree = constructor.build_tree(aln)
Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)

I have an msa file as input and it gives a tree, which looks like:
                                                         _______ Human
                                                 _______|
                                                |       |_______ Chimpanzee
                                         _______|
                                        |       |        _______ Dog
                                 _______|       |_______|
                                |       |               |_______ Cow
                         _______|       |
                        |       |       |_______________________ Elephant
                 _______|       |
                |       |       |_______________________________ Mouse
                |       |
         _______|       |_______________________________________ Platypus
        |       |
        |       |        _______________________________________ Anole_lizard
        |       |_______|
  ______|               |        _______________________________ Chicken
 |      |               |_______|
 |      |                       |_______________________________ Zebra_finch
 |      |
 |      |_______________________________________________________ Xenopus
_|
 |       _______________________________________________________ Zebrafish
 |______|
 |      |_______________________________________________________ Fugu
 |
 |______________________________________________________________ Amphioxus

If you're not wedded to using python to visualise your trees, FigTree provides a great alternative.

Well, this is what hinders me. I am working with python.

Can you provide an example of your input tree? I’d approach this my parsing the treefile itself, truncating the values, writing a new tree, and then plotting that, rather than trying to do it whilst plotting. BioPython’s tree handling is not its strong suit IMO, so for something like this I would probably use ETE3 or Dendropy - those tools may include options for truncating your node/branch values.

I have edited my question for this.

Apologies, try the following,
allclades = list(pars_tree.find_clades(order='level'))
for i, blength in enumerate(allclades):
        printf(blength=%.2f, blength)

"allclades" is a dictionary rather than a simple array.
If this fails I would examine the allclades data structure either using pprint (directly) or using printStruct via the following code,
>>> d = [{'a1':1, 'a2':2, 'a3':3}, [1,2,3], [{'b1':1, 'b2':2}, {'c1':1}], 'd1', 'd2', 'd3']
>>> printStruct(d)

def printStruct(struc, indent=0):
  if isinstance(struc, dict):
    print '  '*indent+'{'
    for key,val in struc.iteritems():
      if isinstance(val, (dict, list, tuple)):
        print '  '*(indent+1) + str(key) + '=> '
        printStruct(val, indent+2)
      else:
        print '  '*(indent+1) + str(key) + '=> ' + str(val)
    print '  '*indent+'}'
  elif isinstance(struc, list):
    print '  '*indent + '['
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ']'
  elif isinstance(struc, tuple):
    print '  '*indent + '('
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ')'
  else: print '  '*indent + str(struc)

Once you know the structure you unravel it accordingly to get to the branch lengths and then issue a printf command to round them to 2 dp.
Goodluck

This still has not helped me. May be an elaborative answer will help.

If I get a moment I will reproduce the analysis, but I'm not your superviser (a very low priority indeed). "list" is the key to the answer in conjunction with printf. Unlocking the object from "list" to get to the branch lengths is the key here.

The universal coding way to do this is the printf command regardless of language (sprintf is used language depending), for Python this is,
blength = 0.22342
printf("Taxa A, blength=%.2f", blength)

Output
Taxa A, 0.22

The following may be useful to your context:
import numpy
allclades = list(pars_tree.find_clades(order='level'))
for blength in enumerate(allclades):
        printf(blength=%.2f, blength)

Goodluck

It did'nt work for me. First there was a syntax error and when I use it this way: allclades = list(pars_tree.find_clades(order='level')) for blength in enumerate(allclades):     print(blength = "%8.2f" % blength) . It says, "TypeError: not all arguments converted during string formatting" Any suggestions?


text:
corrplot issue with correlation plot output
Very conceptual doubt not sure if im doing it right or something terribly wrong so im using this corrplot library to calucate the correlation of my desired data frame after that i rounding it off but when i use my original calculated matrix such as M i get a different plot again when i used the round of matrix i get a different plot in terms of correlation ,in case of M all seems to be positively correlated but when i used the corr i get plot which shows my cell types are both positively and negatively correlated .
Im posting my code below as well as the plot along with the data
library(corrplot)
gendat = read.csv("NEW_RBP/TF_MYELOID_DISEASE.txt",row.names = 1,sep = "\t")

dim(gendat)
head(gendat)
gen <- gendat
head(gen)
M<- cor(gen, use="complete.obs", method="spearman")
corr <- round(cor(M),1)

names(gen)
dim(gen)

mycolors <- rep(NA,23)
names(mycolors) <- names(gen)
mycolors[1:4]   <- 'black' # mpg, cyl, disp, hp
mycolors[5:15]   <- 'red'   # drat, wt, qsec, vs, am
mycolors[16:23] <- 'blue'
mycolors[13:16] <- 'darkgreen'

ord <- corrMatOrder(corr, order="hclust")
newcolours <- mycolors[ord]
newcolours
corrplot(M, tl.col = newcolours, method = "color",is.corr = FALSE,
         tl.srt = 45,diag = TRUE
         ,order = "hclust",hclust.method = c("complete"),addrect = 5,tl.cex = 3,number.cex=1.5)

corrplot(corr, tl.col = newcolours, method = "color",is.corr = TRUE,
         tl.srt = 45,diag = TRUE
         ,order = "hclust",hclust.method = c("complete"),addrect = 3,tl.cex = 3,number.cex=1.5)

I would like to know what is happening in both the cases
data file

corr <- round(cor(M),1) Are you supposed to be calling cor on the result of your previous call to cor? should it just be round(M, 1)?

my goodness that is the mistake im doing.


text:
Are there free cloud computing platforms for biology projects?
I want to implement the analysis found in the paper RNA-Seq of Tumor-Educated Platelets Enables Blood-Based Pan-Cancer, Multiclass, and Molecular Pathway Cancer Diagnostics The project id is 281708. It has 285 samples comprising about 300 Gb of read data. The reads are available on sra.
I want to apply a machine-learning pipeline to find important features in these data but I do not have access to computing resources to finish the project.
Question
Are there free cloud computing platforms that I can use to perform this research?

@Chris_Rands yes it is duplicate.

If you're getting started with a project and need some free computational resources you should look at https://galaxyproject.org/. There are some free servers that you may be able to use. However, Galaxy is designed to have pre-fabricated pipelines and a GUI to avoid having to deal with command-line tools. Therefore any sort of custom analysis will have to go through the galaxy API.
Otherwise, most of the free computational resources that were available several years ago like iPlantCollaborative are no longer available. Running servers costs money, after all, and the funding for such projects often runs out after a few years.
You may consider asking researchers at your local universities to collaborate and provide you with access to their servers.


text:
Database of Cell Volumes by Cell Ontology?
Are there any databases that have aggregated information on cell volumes?  Most useful would be a database keyed by cell ontology (e.g., CL:0000236 for B cell).
The BioNumbers book has a small table, but none of the references appear to link back to systematic catalogues.  In fact, just to put that small table together, they had to reference a dozen or so papers.
The SHOGoiN database is a nice overview of the taxonomy, but the only morphological data available is attached to the images (example B cell), which only includes areas and other histological measurements. Also, they don't seem to provide an API for programmatic access.

I disagree with the "opinion-based" close vote. The question is soliciting users to share any databases that have this information - not to rank them. The examples I gave are in the SO style of sharing what I already tried, neither of which are adequate solutions.

This is a useful question. I searched, and could find nothing matching this.