text:
Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype?
So if I attempt this question using their method I put the data of genotype and phenotype into a data frame in R then use lm() to do linear regression and summary() to get an F statistic with p value 0.54. I'm not connecting what the F statistic means in this case to the answer.
Question 2: Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype? Here, we represent each SNP site as the number of minor alleles on that locus, so 0 and 2 are for major and minor homozygous sites, respectively, and 1 is for the heterozygous sites. We also assume that minor alleles contribute to the phenotype and the effect is additive. In other words, the effect from a minor homozygous site should be twice as large as that from a heterozygous site. You may use any test methods introduced in the chapter. How about permutation tests?

GENOTYPE: 1, 2, 2, 1, 1, 0, 2, 0, 1, 0
PHENOTYPE: 0.53, 0.78, 0.81, -0.23, -0.73, 0.81, 0.27, 2.59, 1.84, 0.03

Answer: Take the linear regression with F-test as an example, the p-value without the permutation test is 0.54 and is 0.55 after 1000 permutation (this number could be different with your answer because of the randomness of the permutation test).
Link to where I'm getting this: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002828#pcbi.1002828-Balding1

What part do you not understand? What did you tried? Did you look at the help pages of the test or any other resource of statistics?

Does the edit make it more clear? Sorry, the first question was poorly written

Yes, much better

Your approach is correct - it's not the only one permitted under their conditions, but it would be the standard way given only that information.
F-test compares a "full model" (predicting the phenotype using genotype and an intercept) versus a "reduced model" (predicting the phenotype using only the intercept, i.e. the mean value of phenotype). The p-value tells whether the full model predicts better, i.e. whether adding the genotype variable significantly improves the prediction. In this single-variable case, it should be approximately equal to the p-value for the beta coefficient of the genotype, if that's more familiar.
Anyway, this isn't specific to genetics, so you can check out other online sources - e.g. https://onlinecourses.science.psu.edu/stat501/node/295 for a longer, but accessible explanation of the F-test, and https://stats.stackexchange.com/ for many answers to similar questions.


text:
IonTorrent variant caller failing, lack of ZM: tags in BAM files
I have a lot of IonTorrent sequencing data in BAM format that I want to call variants on, using the IonTorrent variant caller. However, for some of the BAM files (for about 12 samples out of ~150), the variant caller stops due to what it states is lack of ZM: tags in a few individual reads.
The thing is, most of the reads in the BAM file do not have ZM: tags (when viewed through samtools view), but these do not apparently cause problems.
What settings do I need to change to get these BAM files to process, or do I just filter out the problematic reads and analyse the remainder?

It turns out, after speaking to tech support, that these are old files from a previous version of the software so won't run through variant caller.
The files themselves need reanalysing from the initial FASTQ -> BAM with the newest version.

source: Stack Exchange

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

sort: votes


text:
Why chose CNN for a variant caller
Google released their variant caller DeepVariant which won the highest SNP performance award in the Precision FDA Truth challenge (99.999% accuracy).
From the linked github repo, we see that DeepVariant is a CNN, we provide images of the reads as an input and train the neural network with a stochastic gradient descent algorithm.
To me this is unintuitive. Why use images of reads when you already have the reads data?

https://en.wikipedia.org/wiki/Law_of_the_instrument

Where did you see they need images? The documentation says the input is a bam file, just like any other caller.

@terdon Images (of pileups, of course) are mentioned in the preprint. I think they're also using the Inception CNN, which is heavily used in image processing. I think many people forget that images are just matrices.

@terdon It seems that the first version was on images: We have made a number of improvements to the methodology as well. The biggest change was to move away from RGB-encoded (3-channel) pileup images and instead represent the aligned read data using a multi-channel tensor data layout. and now they use something else. Thank you!

As @DevonRyan has already said, these are not images per se. From a machine's point of view an image is just a bunch of matrices (one per colour channel). I'd say this is just an example of bad/confusing use of terminology.

It's more "images of reads" rather than actual images of reads. In reality, they're feeding in a pileup, which is just a matrix of letters or numbers (you can visualize this as an image, of course). This is more apparent if you look at one of their example python notebooks, CNNs are popular in the image processing field, but at the end of the day images are just large matrices of values so any time you can convert your input into that sort of format then a CNN can be useful (provided you have sufficient training data and some GPUs to throw at things).

I find the answer confusing: are you saying that the CNN is not fed images? Because the way google describes DeepVariant, they are indeed creating pileup image data and putting that through their CNN.

It's not an image, it's an "image". The DeepVariant are a bit ambiguous on that point, though at the end of the day it's a multidimensional matrix of values (the relation of that to an image being obvious).

Why use images of reads when you already have the reads data?

Neural networks take one or multiple N-dimension arrays (n-d arrays) as input; they can't take arbitrary forms of data. To feed reads and their alignment into the network, you have to encode them as n-d arrays. If you ignore "colors" for now, a straightforward choice is a (windowLength, readDepth) 2D-array as you see in an alignment browser – it is an image. The "colors" of this image are often called "channels" or "features" in deep learning (DL). They represent base types and base quality in deepVariant. In the end, you get a (windowLength, readDepth, nFeatures) 3D array as input. This is the most natural encoding to DL researchers.

Why chose CNN for a variant caller

That is because CNN has been repeatedly proved to be the best architecture for such (height, width, nChannels) 3D arrays as input. Read alignments are exactly that.

On a historical note, deepVariant was first built on top of google's now-deprecated DisBelief framework that only takes RGB images as input. They had to squeeze all types of information into three color channels. Now with TensorFlow, an image can have an arbitrary number of color channels. This simplifies implementation and improves accuracy. In addition to this image, deepVariant also takes other information as input.

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


text:
What is the difference between the entity_src_genCategory tag and the entity_src_natCategory tag in pdb-xmls?
I want to read out the taxonomy id of the source organism of a pdb entry using the xml-file of the entry that is provided under http://files.rcsb.org/view/{pdb_accession}-noatom.xml.
I have so far found the taxonomy id in the tag PDBx:entity_src_genCategory (example entry 1UKV) but also in the tag PDBx:entity_src_natCategory (example entry 2BWI).
XML for 1UKV with PDBx:entity_src_genCategory
<PDBx:entity_src_genCategory>
  <PDBx:entity_src_gen entity_id="1" pdbx_src_id="1">
     <PDBx:gene_src_common_name>baker's yeast</PDBx:gene_src_common_name>
     <PDBx:gene_src_genus>Saccharomyces</PDBx:gene_src_genus>
     <PDBx:host_org_genus>Escherichia</PDBx:host_org_genus>
     <PDBx:host_org_species>Escherichia coli</PDBx:host_org_species>
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_gene_src_ncbi_taxonomy_id>4932</PDBx:pdbx_gene_src_ncbi_taxonomy_id>
     <PDBx:pdbx_gene_src_scientific_name>Saccharomyces cerevisiae</PDBx:pdbx_gene_src_scientific_name>
     <PDBx:pdbx_host_org_ncbi_taxonomy_id>469008</PDBx:pdbx_host_org_ncbi_taxonomy_id>
     <PDBx:pdbx_host_org_scientific_name>Escherichia coli BL21(DE3)</PDBx:pdbx_host_org_scientific_name>
     <PDBx:pdbx_host_org_strain>BL21(DE3)</PDBx:pdbx_host_org_strain>
     <PDBx:pdbx_host_org_vector_type>PLASMID</PDBx:pdbx_host_org_vector_type>
     <PDBx:plasmid_name>pET19tev</PDBx:plasmid_name>
  </PDBx:entity_src_gen>

XML for 2BWI with PDBx:entity_src_natCategory
<PDBx:entity_src_natCategory>
  <PDBx:entity_src_nat entity_id="1" pdbx_src_id="1">
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_ncbi_taxonomy_id>223</PDBx:pdbx_ncbi_taxonomy_id>
     <PDBx:pdbx_organism_scientific>ACHROMOBACTER CYCLOCLASTES</PDBx:pdbx_organism_scientific>
  </PDBx:entity_src_nat>

What is the difference between them? Are there other tags that contain information about the source organism of the entry? Can they occur together in one file? If so, is one of them more valid?

src_nat = isolated directly from a natural tissue
src_gen = the source was genetically manipulated
The files that you use (PDBML/XML) are derived from the primary PDB data format which is PDBx/mmCIF.
The names in XML seem to correspond to mmCIF names, so you can read documentation on http://mmcif.wwpdb.org/, and here are direct links:

category entity_src_nat
category entity_src_gen

Can they occur together in one file? I suppose they can for complexes: one "entity" could be natural and another genetically manipulated.

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


text:
XCMS isn't evaluating the proper peak due to retention time drift
I'd like to do a targeted screening analysis for trace contaminates.  However, due to the nature of the sample preparations, the retention times I'm getting drift a lot from run to run and XCMS isn't properly identifying my standards and thus affecting my ability to evaluate my samples.
Is there some trick to getting XCMS to identify the proper retention time that I'm missing?
Is there some other software I should be using? (We have several different instruments from different vendors so it's better if it's not tied to the particular MS)

I'm not an expert of mass spectrometry, but perhaps showing the figures or some data might provide an expert with some clues to help you.

As Llopis said, some more details should help to attract an answer. Particularly, which version of XCMS and showing XIC's from some of your standards to show the run-run variation.

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


text:
A good tool for gene locus visualization
Does anyone know of a good tool for visualizing gene loci? Basically, I'd like to show how a group of genes look in different strains of a bacteria. I have annotation of course, I just need a tool to plug it into and get an image of the locus to compare.
I'm talking something basic like this:

I ended up using Easyfig, which is made for the task and requires little tinkering around, once fed the genbank files.
(Granted, being so simple, it's hard to make the image more advanced, but it serves my purposes for now).

Here's an example of my image.

I think in R, Gviz is the best choice. I have used it for human sequences but it will probably also be suitable for Prokaryotes.
Another option is MS powerpoint, if you only need the arrows...

yah +1 gviz is great, i've used it for bacteria

I wanted to avoid putting something together manually with powerpoint, even though I suspect that's how my supervisor does it...
I'll take a look at Gviz though!

GFF2PS was developed for exactly this sort of thing. It takes a GFF file as input, and will produce a pretty PS output showing the annotations. It was used to create the poster for Celera's human genome draft [1].
Here's a very simple example of what it can do:

And a more complex one from the mosquito genome paper [2] (click here for a full size, 1.5M gzipped PostScript file):

If you are happy with the defaults, it is very easy to use, but otherwise it is a bit of a pain to configure. It is absurdly powerful though.

Venter et al. Science 291(5507):1304-1351 (2001)
Holt et al. Science 298(5591):129-149 (2002)

gggenes is a ggplot2 based solution and looks very appealing. Plus, it has a nice introduction with examples.
library(ggplot2)
library(gggenes)
ggplot(example_genes, aes(xmin = start, xmax = end, y = molecule, fill = gene)) +
  geom_gene_arrow() +
  facet_wrap(~ molecule, scales = "free", ncol = 1) +
  scale_fill_brewer(palette = "Set3") +
  theme_genes()


text:
Error while Importing fastq sequences to jMHC software
I am trying to analyze my .fastq files (MiSeq exon reads) with the jMHC software, but output files I have obtained have no results.
I think I might be doing something wrong during the first step: "connect to database".
Because on the "import/extracting" step, I open all my .fastq files, click "Start", and it shows this message:

"Success! Processed: 139, Imported: 0".

So I assume it is reading the files OK, but not importing them to the "database"? Could it be the extension of my files? The software says it operates on FASTA files. Are .fastq files not compatible?
Anybody that has successfully run this program could give me a hand?

Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.

Hello, thanks for the answer. It would be useful if you would provide link to the tool you suggest.


text:
Volcano plot in R
This question has also been asked on biostars
How can I reproduce this volcano plot?

I'm only able to do the traditional one, I'm kind knew too these field.

This plot is clearly done using core R functions. There are smoother alternatives how to make a pretty volcano plot (like ggplot with example here), but if you really wish to, here is my attempt to reproduce it :

I obviously had to generate data since I do not have the expression data from the figure, but the procedure will be about the same with the real data. Here is the R script :
### define interesting genes
interesting_genes <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
fc_interesting_genes <- c(-1.8, -1.3, -0.9, -0.85)
pval_interesting_genes <- log10(c(0.0005, 0.01, 1e-4, 4e-05))

### define colours (grey, blue, red)
palette <- c(rgb(0, 0, 0, max = 255, alpha = 65),
    rgb(0, 0, 255, max = 255, alpha = 125),
    rgb(255, 0, 0, max = 255, alpha = 125))

### generate some expression like data, fold change using normal distr, and tranformed uniform for pvalues
fold_changes <- c(rnorm(1000, 0, 2), fc_interesting_genes)
pvalues <- c(log10(runif(1000) / 3^(1 + abs(fold_changes[1:1000]))), pval_interesting_genes)

# find deferentially expressed genes to colour
left_biased <- fold_changes < -1
right_biased <- fold_changes > 1

point_type = 20
point_size = 1.5

### plot corpus
plot(fold_changes, pvalues,
     pch = point_type, cex = point_size, col = palette[1],
     xlab = "Fold Change (log2)", ylab = "P value",
     xlim = c(-6,6), ylim = c(-5, 0), yaxt="n")
aty <- axTicks(2)
labels <- sapply(aty,function(i)
            as.expression(bquote(10^ .(i)))
          )
labels[6] <- 1
axis(2,at=aty,labels=labels)

### plot deferentially expressed genes in colour + genes on interest with black
points(fold_changes[left_biased], pvalues[left_biased],
       pch = point_type, cex = point_size, col = palette[2])
points(fold_changes[right_biased], pvalues[right_biased],
       pch = point_type, cex = point_size, col = palette[3])
points(fc_interesting_genes, pval_interesting_genes, pch = point_type, cex = point_size)

### add text
labelx <- c(-4, -4, 4, 4)
labely <- c(-1.5, -0.5, -0.5, -1.5)
labels <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
text(labelx, labely, labels)
text(-4, -4, "down\n194 genes")
text(4, -4, "up\n156 genes")

### add lines between text and data points
for(i in 1:4){
    lines(c(labelx[i], fc_interesting_genes[i]), c(labely[i], pval_interesting_genes[i]))
}

Without a minimal reproducible example (MWE), I can't reproduce the plot but I would suggest using existing volcanoplot functions such as the limma package on Bioconductor. Here is a typical workflow for a differential expression analysis that produces a violinplot:
#install package
source("https://bioconductor.org/biocLite.R")
biocLite("limma")

#load package
library("limma")

#define phenotype/experimental conditions for samples
group = rep(c("Treatment","Control"),
c(number_of_treated_samples,number_of_controls))
design = model.matrix(~group);
colnames(design) = c("Control","TreatmentvsControl")

#fit empirical Bayes model
fit <- lmFit(NormalisedExpression, design)
fit <- eBayes(fit)

#summarise significant results
tt <- topTable(fit, coef="TreatmentvsControl",
                adjust="BH",n=nrow(NormalisedExpression))

#plot results
volcanoplot(fit, coef="TreatmentvsControl")

This plot can be customised in a similar manner to base R plots by passing the relevant arguments as shown in the limma documentation. This is the recommended plot format that readers in the field will be familiar with. While I do not recommend deviating from this standard, it is possible to reverse axes in R and the source code for the plot function can be retrieved from R by entering volcanoplot (without any arguments) if you wish to modify it.

Easy to visualize volcano plot using Python. Check this tutorial for inverted volcano plot


text:
How to read gene regulatory network edge list files?
The following is an excerpt from an edge list file from the Gene Regulatory Database, YeastNet v3
YML100W YMR261C 5.73968189451551
YDR074W YML100W 5.73143623470086
YDR074W YMR261C 5.72385204041689
YML106W YMR271C 5.64698433340441
YGR240C YMR205C 5.64063955646679
YDL131W YDL182W 5.63907349642605
YNL037C YOR136W 5.62782704009976
YHR208W YJR148W 5.61712512462848
YAL062W YOR375C 5.59956647235823
YNL104C YOR108W 5.59046546303091
YGR132C YGR231C 5.58996626839351
YBR140C YOL081W 5.58738025271142
YDR453C YML028W 5.58281825879471
YNL096C YOR096W 5.58060076040965
YGR180C YJL026W 5.57785446902287
YGR056W YLR357W 5.57582611528462
YBR009C YNL030W 5.57261292043205
YMR169C YMR170C 5.55866364259678
YAR071W YBR093C 5.55668852247789
YER081W YIL074C 5.55252000887001

What do these rows mean?

It also depends which file you have downloaded. Apparently the evidence can be based on various inputs such as co-citation in PubMed, or co-occurrence of domains... The details are listed on their webpage.

You can read about it on the website, the last column is the Log likelihood scores (LLS) between the genes in the first two columns. An LLS of 0 means not better than random link.
On the website they show papers with fuller details about the statistics.


text:
Seurat Merged objects tSNE - How to paint on original IDs?
I am working with single-cell RNA-seq data, using the R package "Seurat" to cluster and visual data-points.
I had two single cell datasets from which I generated two Seurat objects. I then combined the two using MergeSeurat. I did differential gene expression analysis, performed clustering, and ran a tSNE plot.
Now what I want to do is create that same tSNE plot, but color each cell according to the original object from which it came. How do I do this?

@MrFlick Don’t you agree that it’s a fairly clear question for domain experts (which, according to your profile, you are)? I’ve flagged it to be moved to bioinformatics.SE.

Below a few lines of code that accompany BC Wang's answer.
After using MergeSeurat the sample name will be added to meta data under orig.ident. this can then be used to color the tSNE either using group.by or pt.shape. The former will show colors for each sample, the latter will color each cluster and give sample id another shape.
path1 <- file.path("path to 10x sample 1")
path2 <- file.path("path to 10x sample 2")

S1 <- Read10X(data.dir = path1)
S2 <- Read10X(data.dir = path2)
S1 <- CreateSeuratObject(raw.data = S1, min.cells = 3, project = "S1")
S2 <- CreateSeuratObject(raw.data = S2, min.cells = 3, project = "S2")

s.combined <- MergeSeurat(object1 = S1, object2 = S2, add.cell.id1 = "S1",
            add.cell.id2 = "S2", project = "2samples")

All data analysis goes here:
RunTSNE etc...
Plot the data
TSNEPlot(s.combined, group.by ="orig.ident")
TSNEPlot(s.combined, pt.shape="orig.ident")

If I understand your question correctly, you are trying to present the tSNE plot coloring and labeling the cells representations based on the original sample ID, i.e. which Seurat object it was from.
I think they offered a vignette on that. Please check this link. https://satijalab.org/seurat/merge_vignette.html
Basically, you can add the sample name for the cells in the argument option in MergeSeurat function. Please check the documentation on MergeSeurat function as well.
When I did this, I have a lot of samples to be merged, so I simply saved the sample info for every cell in each sample to different objects. You can also use AddMetadata to do it as well.
Hope that helps.

It is impossible to create the same tSNE plot without knowing which seed you used. tSNE plots are not like PCAs you don't get the same results when using the same data.
But if you want another tSNE plot with a color differentiating from which object each sample comes you'll need a column with these information and color the plot by it.

Seurat allows you to visualize tSNE plots with gene expression data painted on top using the FeaturePlot function. I assume there must be a way to do the same thing with meta data such as original IDs. I just can't find anything online that shows how to do this.


text:
Seurat for clustering bulk RNA-seq?
Is it ever ok to use Seurat for clustering bulk samples?
I am looking at FPKM data from ~750 bulk RNA-seq samples generated using Cufflinks. As suggested for FPKM data, I manually input log transformed data to the @data slot [cd138_bm@data <- log(cd138_bm@raw.data + 1)] and skip the NormalizeData() function. I then use functions FindVariableGenes, ScaleData, RunPCA, FindClusters, RunTSNE, FindAllMarkers in their usual ways to find clusters & cluster markers. My clustering results are quite reasonable, and reflect published work clustering similar samples.
What are the potential pitfalls of using these Seurat functions on bulk data? In FindAllMarkers, would you recommend I use the "negbinom" test? (currently using wilcox) Any other arguments you would recommend changing from the default?

Good question. I removed scrnaseq tag since you have bulk rna sequencing, if I misunderstood somehting, feel free to [edit] it back.

I also want to try this bulk analysis with seurat and I am interested in your methods.
Did you publish this analysis already ?
If it's fine, I would like to ask you get the papers you wrote !

I didn't publish, but I did use seurat on bulk and it worked great.

It's a good question. I processed the sc-RNAseq data and generated some new markers based on gene expression value. In the subsequent analysis, I want to find differential markers in multiple clusters of my dataset based on seurat package.

I'm not sure Seurat is the best tool for this as it was developed for single cell RNA seq data and there are a few intricacies of that type of data that are very different from bulk RNA seq.
For bulk there are really good packages available and corresponding workflows, e.g. limma, edgeR and DESeq2.
The main problems with Seurat for bulk RNA-seq:

Seurat expects counts as input - FPKM are not counts nor are they log counts or log norm counts. It's a specific normalisation method that takes into account gene length and library size and breaks the link between gene counts and variance.
Because of the distribution of reads from a single cell that you get from the sequencer Seurat scales all unique reads (i.e. UMIs) to 10,000 per cell and usually regresses out the nUMI count. Again, becuase you've got FPKMs this makes no sense.
The other functions you can use on bulk data, i.e. tSNE, finding clusters and markers (which is essentially DGE analysis). The one thing to note for the latter is that again because you've got FPKMs they don't follow the negative binomial distribution. Your best bet is limma trend but I'm not sure it's included in Seurat.

You don't need to use Seurat's normalization. You can provide normalized counts or FPKMs and skip the normalization step. See this earlier discussion: https://bioinformatics.stackexchange.com/questions/5115/seurat-with-normalized-count-matrix

That's true but if you look at the original post, OP uses ScaleData, which scales already scaled FPKM values. Even though feeding bulk RNA seq FPKM values to a Seurat@data slot is essentially doable, you are still relying on the downstream processing steps that were designed for log-normalised count data to be meaningful for FPKM values and it's certainly not the case for the FindAllMarkers function.

There are a few issues on the Seurat GitHub that discuss importing TPM/FPKM values. The suggestions are always to just skip normalization and proceed normally. I haven't seen any concerns raised about FindAllMarkers. I don't mean to imply those concerns are not legitimate, but it's interesting they are not raised.

FYI - heard from the author of Seurat, see this github issue
His suggestion to use DESeq won't work in my case, because DESeq requires counts data which I don't have. But it is the green light to generally use the package on bulk data. I am using test.use = "wilcox" in FindMarkers() because it's non-parametric...

I would actually suggest limma-trend (not voom). Compared to a wilcox you'll gain a bit of statistical power due to the empirical bayes moderated t-test.


text:
Separation of mixed plasmid DNA sequences post whole-plasmid sequencing
Imagine a DNA sample containing a mixture of different intact plasmids. These samples are sequenced using either MiSeq or HiSeq sequencing. Would it possible to assemble these plasmids post-sequencing as they would have been when sequenced individually?

Would you know the sequences? I mean, will you have reference sequences for each plasmid which you can use for assembling? How similar are the plasmid sequences? Will you have multiple ~100nt sequences that align perfectly to many different plasmids? Please [edit] your question and add more details.

Don't forget a template assembly, to confirm your results once you are clear what plasmids you are dealing with .. e.g. bowtie2, any paper on bacterial metagenomics will help you. Read length vs. sequence similarity is important

You would expect to have high coverage, given the plasmids are short, so de novo assembly would be likely very easy.  Given that each plasmid is present in different multiples, you would expect different coverage on each plasmid, so it might be best to approach it as a metagenome-type or transcriptome assembly, rather than a classic genome assembly.
Alternatively, Google search just told me there is a tool called plasmidSPADES designed for just your purpose.  SPADES itself is a reliable genome assembler, so I would imagine plasmidSPADES to be a good choice.  It does assume that there is also some whole bacterial genomes in the sequencing mix though.
As @terdon commented, it does depend on whether you have many identical copies of each plasmid, or many similar but not identical plasmids in the mix.

Yes I had the same ideas. I am very acquainted with SPADES,although never tried de novo assembly on a mixture of plasmids.

It looks as though plasmidSPADES expects to find some 'long' chromosomes in the mix, so it might be important to throw in some whole genome DNA into your mix of plasmids, to train the assembler.  You could do this after the experiment, by appending a fastq file of E.coli genomic DNA to your fastq file of plasmid sequences...

Coincidentally, the DNA sample is "contaminated" genomic DNA since it was extracted from an environmental sample. That should provide the necessary SPADES training?

If it's from an environmental sample might be better to try a metagenome assembler, maybe metaspades.  Alternatively, use kmer spectrum to try to pull the dataset apart based on abundance, and run the assembler on cleaned-up subsets of the reads

Let's take a step back and consider the "perfect" output for a de novo assembly algorithm. Ideally, you would like to see one complete sequence for molecule (chromosome, plasmid, etc.). In reality, this is difficult to achieve due to a couple factors.

By random chance some regions may have low coverage, meaning that there are an insufficient number of reads spanning the region to reconstruct it de novo.
A bigger problem is that some portions of the genome are repetitive and occur at many different locations, sometimes even on different chromosomes.

So in practice, it's very common that de novo assemblies can only recover partial fragments of each molecule, and sometimes fragments are erroneously fused due to shared repetitive content.
Back to your question: it's certainly possible to assemble distinct plasmids that were sequenced together without any special sample prep. How well this will work in practice will depend on the amount of coverage/redundancy in the sequenced reads and the amount of DNA shared between the plasmids. If they don't share much in common, and you sample to a sufficient depth of coverage, you shouldn't have any problem assembling them.

I think your comments are very informative so thank you so much!

@RoelofCoertze I'm glad it was helpful!


text:
Using Seurat to compare mutant vs.wt
I am interested in using Seurat to compare wild type vs Mutant. I don't know how to use the package. How can I test whether mutant mice, that have deleted gene, cluster together?

This question is too vague and open-ended for anyone to give you specific help, right now. Can you provide a bit more detail? Have you run Seurat on your WT samples and have well-defined clusters that you want to compare against in your mutant samples?
If you can show some figures or briefly describe what you've done so far, or specifically what you're trying to achieve, it will be much easier for someone to provide guidance.

Single-cell analysis to compare samples is a long a difficult process. There is very good documentation for 10x Genomics cellranger, the DropSeq Pipeline and the Seurat R package. These tools all have GitHub repositories and the authors are very responsive if you encounter issues.
Depending on the technology used to generate the data, you'll need to use either cellranger or DropSeq to process the FASTQ files. These are designed to account for the the different experimental design. Still, it is important to ensure that your samples have been processed in the same way to reduce batch effects.
This question is to open-ended to cover all of the details needed to perform the analysis but here is some of the main things to consider and the general process to perform such as analysis. You will need competence using the command-line (shell) and programming in R (or Python) to perform single-cell analysis. I recommend performing this analysis on a remote server as some steps are memory intensive.
The general overview to compare samples is:

Demultiplex the reads for different indicies. You need to do this even if you've only put one sample per lane (to remove the index from the reads). cellranger mkfastq or Illumina's bcl2fastq will do this.
Obtain a reference genome (FASTA) and gene annotation (GTF) for the species you are working with. You can prepare a reference transcriptome with cellranger mkgtf and cellranger mkref.
Run cellranger count or the DropSeq pipeline on each sample separately. These will both perform STAR (splice-aware) alignment of paired-end RNA-Seq reads and count UMIs for each cell-barcode and gene. Cells will be identified and filtered by an automatic UMI threshold (this can be changed by forcing with expected number of cells). This will return a gene-barcode matrix.
You need to aggregate results of cellranger runs for different samples with cellranger aggr. This will perform downsampling by default to normalise the number of reads between samples. It will also perform graph-based (Louvian) clustering on the combined samples and return HTML and Cloupe summary data to explore, in addition to a combined gene-barcode matrix for all samples.
Single-cell RNA-Seq experiments are subject to batch effects. The CCA method from Seurat and the MNN method from scran are both available in the respective R packages to account for this (with different approaches). Batch effect correction needs to be performed before downstream analysis to ensure comparisons are valid.

This methods require an overlap between the samples to estimate technical errors and account for them. You need to bear in mind that neither of these methods works unless:

some of the cells detected in each sample are expected to be the same cell type
biological replicates of each sample have been performed

Once you have performed this correction, then standard Seurat workflows such as those shown in the tutorials can be performed which will identify clusters. You can use these to identify clusters specific to each genotype. Bear in mind that some overlap should be expected and if there isn't one, batch effect correction will over-correct for this.

Hi Tom, Thank you for your advice. You have mentioned that we can use a remote server. Yes I have tried that. But I don't know how to save our data in a remote server. Can you tell me how to? Thank you so much!

CCA corrects based on shared variability. ie. if you expect to find a unique new cell type or cell state for your mutant condition than after CCA that cell type will more often than not be blended in with the other clusters. You won't detect it. If this is what you expect MNN might be better. However if you expect changes in gene expression between the same cluster for wt and mutant CCA might work just fine. Those limitations have been addressed very recently by many groups.

MNN has similar limitations. It’s stated as such in the original publication of this method.

You just need to run cellranger, it will write output files to disk. You save data to disk from R on a remote server. This is really a separate question and I think you’ll find answers for how to do it elsewhere. You need to ask specific questions in SE to get the best answers.

There is no single ended answer to your question as everything will depend on your data. For example: Is this unsupervised data where you have to browse to 100 cell types for each sample or did you FACS sort on specific markers and are you looking for specific differences between 2 specific cell types etc...
Below you can find 4 recent packages that assist in your question. I strongly suggest you read the preprint of each, play around with your data and then follow the suggestion in the comments and ask more specific questions about how to proceed correctly.
If you want to use Seurat you can start here:
Seurat 3.0 pre-release to integrate your data sets:
https://satijalab.org/seurat/pancreas_integration_label_transfer.html
Seurat 2 tutorial on CCA to compare different conditions
you can use these methods on your integrated data sets.
https://satijalab.org/seurat/immune_alignment.html
LIGER Package and preprint can be found here
https://macoskolab.github.io/liger/
HARMONY
https://www.biorxiv.org/content/early/2018/11/05/461954
Conos
https://github.com/hms-dbmi/conos


text:
How to generate Genbank format out of paired end fastq data?
I've unmapped cleaned paired end sequence data in fastq format of a bacterial genome. I want to get a sequence data in Genbank format in the end. What are the exact steps that I have to follow in galaxy webserver and if in some point necessary, also using linux based apps? I would appreciate a lot of details.

The steps needed to assemble a genome in Galaxy are layed out in the Assembly training material from the Galaxy Training network. The software for this should be available on the main galaxy servers (e.g., usegalaxy.eu). Since you have a prokaryote, I suggest you follow the "Unicycler Assembly" workflow, which is also available on usegalaxy.eu.

I think the reason for your confusion is vocabulary.
The keywords you search for are genome assembly (see @Devon Ryan's answer) a process that is trying to stick together reads to recreate the sequenced genomic sequence. The assembled sequence needs to be annotated if you are interested where are the genes or what they could possibly do. If you will search for tutorials using these keywords you will find plenty.
I am not sure why you need it, but Genbank format contains the genome (the product of genome assembly), metadata (who is the submitter, what species it is, etc) and genome annotation (what is in the region of genome).


text:
Get a list of PDB ids from a list of sequence protein number?
I am working on a project about phosphorylation, and I am creating a database of some proteins. I have a list of protein accession ids for which I need to get the corresponding PDB file for each one to map it to PSM matrix. I want to get the PDB file using Python or BioPython tools. Any suggestions will be appreciated.

Thanks for your reply. I have a list of protein ids, and I want to retrieve PDB code for each one! I made the protein ids database with the phosphorylation cite for each position. However, I need PDB code because I want to map it to PSM matrix. I have no idea how to retrieve PDB code! Thanks again

Please include a few IDs and the expected outcome, i.e. ID: P26378 should map to PDB: 1FXL, etc. If you just need to do it a few times, you could use the UniProt web interface, just upload a bunch of identifiers and get the corresponding PDB IDs. https://www.uniprot.org/uploadlists/

Thanks man, I really appreciate your reply. I am currently using this website to map protein accession ids to  PDB. However, I want to do it by the code because I have a big list of protein ids. I want to get PDB by code so I can map it to PSSM to get the structure of the proteins. I want to do it in the one piece of code,not separately! Again Thanks for your kind and helpful reply!

A typical Uniprot URL would look like the example here:
https://www.uniprot.org/uniprot/?query=ID:P26378&format=tab&columns=id%2Cdatabase(PDB)
The response will be:

Entry Cross-reference (PDB)
P26378    1FXL;1G2E;

We want our query to include the protein ID, the output format should be tab separated and the only contain the columns id and database (PDB).
To get the PDB IDs you would need to split the response, take the last line, split this line after each tab and split it again after each semi-colon.
Using Python's Requests library the whole code would be.
import requests
import time

uniprot_ids = ['P26378', 'O35433', 'Q02910']
url = 'https://www.uniprot.org/uniprot/'

protein_to_pdb = {}
for protein in uniprot_ids:
    params = {
        'format': 'tab',
        'query': 'ID:{}'.format(protein),
        'columns': 'id,database(PDB)'
    }
    contact = ""  # Please set your email address here.
    headers = {'User-Agent': 'Python {}'.format(contact)}
    r = requests.get(url, params=params, headers=headers)

    protein_to_pdb[protein] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')
    protein_to_pdb[protein].pop(-1)
    time.sleep(1)  # be respectful and don't overwhelm the server with requests

print(protein_to_pdb)

To accelerate things, you could store your protein_to_pdb dict in a file and only retrieve the new ones
Checking if the returned ID matches your input ID is always a good idea

Thanks man, I really appreciate your answer and your time with me. Please I want to read a file which has a list of protein ids as df ,and retrieve all PDB codes for each one. I am trying to edit on your code by using :"df.iloc[i,1] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')" ,but there is an error which is "list index out of range". Please can you  help me with that?

I did it by using :"df.iloc[i,1] = r" instead of using a list,and I get a file with all protein ids and PDB codes. Is that correct? Or I am missing something! Thank you so much for your answer; it is really helping me a lot. I really appreciate it .

Try a more simple example, just copy a few ids in the list and see if it works. Without seeing your dataframe it is hard to guess where the problem occurs.

You could use mapping between UniProt and PDB entries from the SIFTS project.
The uniprot_pdb.csv file there has UniProt accession names followed by a semicolon-separated list of PDB four letter codes.


text:
how do i change a variable with iterator under a for loop in python
This question comes up in my bioinfo class.
I have a csv file containing gene names and the cluster they belong to. I want to run GESA on each individual cluster. So I would like to extract the list of gene names corresponding to each cluster. So I wrote code like this
import pandas as pd
gene_list = pd.read_csv('marker_genes.csv')
for i in range(12):
    "cluster_"+str(i) = gene_list.gene[gene_list.cluster==i]

I know this raises syntax error, so what is the correct way of doing it.
In the end, I would want cluster_0 to cluster_12 each being a list of gene names.

You're looking for the groupby() function https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html

Groupby is probably right, not always easy to get the object to spit out the values. Difficult without knowing the table format

Reshaping via melt is another possible way in or pivoting might work. There's probably lots of solutions. The pd.read_csv is a neat part of pandas, however you rarely use loops, only when you are stuck and its a different syntax.

cluster_0 to cluster_12 are just strings. they cannot be the names of arrays or matrices. Strings can only be used in print statements and string processing statements.
The correct way to define an array or matrix is as follows: call the arrays as: cluster[0] to cluster[12]
then, assign:
  for i in range(12):
        cluster[i] = gene_list.gene[gene_list.cluster==i]

then whenever you need to print the names of these arrays, you use:  "cluster_"+str(i)

Okay .... yeah, this is better. There is an issue with the "cluster" dictionary (??), so I still think it will throw an error, but it is a definite improvement. The pandas command that is the formal solution is "iterrows" for this style of loop. – its a fun question.


text:
Unmapping Alternative reads in BAM file
I want to do some HLA typing, most of the tools require the bam file is aligned to primary genome without alternative read handling. From the International Genome Sample Resource project, I can get the high coverage file, but they are all aligned to primary genome with hla, decoy and alternative reads. I want to make the alternative reads into unmapped reads.
I think about unmapping the whole file to fastq, and then realigning according to different assembly, but it's time-consuming and also costs a lot of storage, so I don't want to follow that approach.
I also tried this script, but so far it's not working (I have all requirement packages installed).
Any help, tool or insight for this problem is warmly welcome.
Thank you

hmmm are you familiar with samtools? you can convert unmapped reads on the fly (on the console) and then remap it and output a bam file

the only thing i am not sure about where do your HLA reads go in the genome, if they are multimapping or you know that they map to a specific location, you can also do something similar to above

thank you for your comment. I'm not so familiar with samtools, actually, any answer for me is valuable. For the HLA allele, I use package like optitype and kourami, they can handle it if I can give them the correct input.


text:
CollectHsMetrics base coverage output has overlapping targets
I am trying to get base coverage information, and am using the --PER_BASE_COVERAGE output from Picard tools to get a text file of base level coverage. However, some of my targets are being grouped together. In the targets file, these are separated by one base as such:
12      12258595        12258695        +       TARGET1
12      12258696        12258795        +       TARGET2
12      12258796        12259791        +       TARGET3

However, when I look at my covg file, I see the following:
chr    pos    target   covg
12  12258595        TARGET1|TARGET2|TARGET3 510

Why would HSmetrics group these targets together, and how can I separate them?

Question is also asked on biostars.

Hello @Blaze9, I don't understand your question. Would you maybe explain bit what are the targets? Are they read groups? (disclaimer, I am not a Picard tools user); I usually use samtools for this.

Hi @KamilSJaron, the targets are either tiling regions or exons/introns. I can checkout samtools to see if it provides better results. The per_target_coverage output does not combine the targets which is the weird part. Thanks!

Ah, I see. Maybe you want to [edit] your question and add as many details as you can (such those you wrote in the comment) :-).


text:
Where can I find .cdf files for Affymetrix ChIP-chip data analysis?
I am analysing microarray ChIP data generated several years ago. I want to use the affy package for R. However, in order to load the CEL files and proceed with the analysis, I need the pombe annotation file which comes under the name of Sp20b_M_v04.cdf. I have looked everywhere for this file but I am unable to locate it. Is anyone in possession of this file or know where I can download it from?

Welcome to Bioinformatics.StackExchange @Min! Could you let us know where the instructions indicate that you need the Sp20b_M_v04.cdf file? Is this from a colleague's lab notebook, or from the affy package's documentation, or somewhere else? This info will be helpful in troublehshooting.

Hi Min, completely agree with @DanielStandage. You need to [edit] your question and add lot more details, otherwise, it's really hard to help.

@KamilSJaron I converted your answer to a comment since it wasn't really providing a full answer to the question.

@terdon I was about to delete it because it was the CEL file I linked, not cdf. ^^

@min I found the corresponding ncbi geo page. it seems that the file you are referring to isn't uploaded. You must contact authors I believe.

Crossposted with biostars?

It looks like it's in GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10187
There is apparently a zipped copy linked at the bottom in the supplementary files.
When people upload a study based on an array, there is a corresponding 'platform' record where the array details should be deposited.
Edit: you will probably need to gunzip the file first, to get it into CDF format.


text:
Retrieving Metabolic networks from WIT Database
How can I get the dataset for all the 43 metabolic networks for this paper:
https://www.nature.com/articles/35036627
It looks like the WIT database does not exist anymore.

It seems like there are some files related to the databases at the ftp site: http://ftp.mcs.anl.gov/pub/WIT2/. The related files were last updated on 2002.
I found it via searching at google: site:*.mcs.anl.gov/ WIT.
Which reports also a page about the database. There seems to be a project to develop WIT3, because there is a mailing list for developers, maybe you can ask them more details if you have problems with the database.


text:
Unable to install bedtools on windows 10 ubuntu
I am trying to install bedtools on windows 10, but I get an error I don't understand. How can I fix it?
Building BEDTools:
=========================================================
DETECTED_VERSION = v2.27.1
CURRENT_VERSION  =
Updating version file.
 * Creating BamTools API
- Building in src/utils/FileRecordTools
  * compiling FileRecordMgr.cpp
make[1]: g++: Command not found
Makefile:36: recipe for target '../../../obj//FileRecordMgr.o' failed
make[1]: *** [../../../obj//FileRecordMgr.o] Error 127
Makefile:186: recipe for target 'src/utils/FileRecordTools' failed
make: *** [src/utils/FileRecordTools] Error 2

Hi there, I suppose you use the bash terminal available at windows 10, is that correct? The error says that you don't have a compiler for c++ called g++. You can try install it by sudo apt-get install g++. Then run make again.

Turns out that something went wrong the first time I tried to install g++, but I didn't see the error message.

My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues.
The linux version of conda installs just fine in WSL and bioconda hae pretty much any bioinformatics package you might need, already compiled.

make1: g++: Command not found

You do not have g++ installed on your machine.  This is used to compile the bedtools into an executable file.
Here is a tutorial on installing g++

If you have admin privileges, the easiest way to install g++ is sudo apt-get install -y build-essential.

@DanielStandage does that work on Windows? I assume the OP is using the newfangled Ubuntu Subsystem for Windows thing, but I don't know if that can install arbitrary packages like that.

Oops, I missed the "Windows" part and only saw Ubuntu.  It looks like apt-get is supported though. https://www.howtogeek.com/261449/how-to-install-linux-software-in-windows-10s-ubuntu-bash-shell/

I had tried that, but another error popped up without it realizing. However, the following error popped up after entering make, and I can't tell if this is critical for the functionality of bedtools:
`- Building main bedtools binary.
done.

Creating executables for old CLI.

/bin/bash: python: command not found
Makefile:145: recipe for target 'all' failed
make: *** [all] Error 127`

Please add the new issue to original post.  This is saying that you don't have python installed.  I think that if you don't want to use the legacy CLI (command line interface) you should be all set.  You can try running bedtools and see if you get a usage prompt.


text:
Will using smaller kmers help get larger contigs? If not, then what?
I've been using MEGAHIT to assemble metagenomes, with particular focus on specific genomic areas.
Sometimes all I get is gene fragments or pathway fragments (eg. if I know that genes A, B, C, D and E should be together, I only get A, B, C in one contig, and maybe D and E in another). That is understandable when the sequencing depth is low, and I'm assume that it's the best I can get.
However. In two of my metagenomes, the genes of interest seem to be abundant (depth of 100-300+ listed for MEGAHIT contigs). And I still have the same problem. I assume that these 'breaks' between contigs are due to natural sequence variation which breaks the assembly process. (while sequencing errors are also possible, I've used bbduk to clean my data prior to assembly)
I'd like to have larger contigs because I'm curious about the gene order. Things I've tried (did not work): visualizing assembly graphs with Bandage around Blast hits, starting with a smaller k-min (21 as opposed to 27).
So I guess I'm wondering if using even smaller k-mers might make a difference. Or if someone has another suggestion, I'm happy to listen! Thank you.
Megahit setting used:
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 1500 --presets meta-sensitive -o output

Meta-large starts from minimum kmer size 27; meta-sensitive from 21.

Typically the tradeoff is smaller kmers getting you a more connected assembly graph, while longer kmers can span repeats more and perhaps connect more contigs. So I think the logical thing might be to try and increase kmer size.

@NatWH Thank you, I can see how that would work. Considering that a MEGAHIT run goes all up to k-mer size 127 (maximum), and the results that I get are as such, shall I assume that nothing can be done just by changing k-mer size? (or even in general)

Does MEGAHIT build the results of the graph from larger kmer sizes on top of those from smaller in some way, akin to SPAdes? If so, then the method should already have accounted for this. If not, then perhaps you can try and manually increase kmer size.

You can try s-aligner (free for 15 days) which usually gets quite larger contigs and quite larger NG50 for metagenomic data containing viruses (also phages). Indeed metaSPAdes also gets larger contigs, despite shorter than s-aligner.

Ref. s-aligner: a greedy algorithm for non-greedy de novo genome assembly

Do I understand right that you are the developer of this commercial product? Please, disclose your conflict of interest.


text:
Extracting expression data in r from GEO and got a S4 object
I want to open the expression data GSE9838 dataset from GEO using GEOquery R package.
gset <- getGEO(filename="GSE9838_family.soft.gz")

The result of this code said it's S4 object of class GSE.
How can I get the expression data of GSE9838?

Welcome to the site. Have you looked for the tutorial of GEOquery? You can see them with vignette(package="GEOquery") . Also have you looked around on this site and in support.bioconductor.org for how have other people done this?

Possible duplicate of Extracting expression data from GSE dataset downloaded from GEO

Not sure if that method works in this case

@Oka Why shouldn't it work? Have Qiran tested it? At least mentioning that it didn't work it would help the people trying to answer the question (and future readers)

@llrs I´ve checked if it would, and it didn´t (at least on the first trial). The method worked fine with GSE used in the example, but didn´t seem to work in this case for some reason..

Ok, how it didn't work, could you post the code so I could see what is the error message and the options you used?

datExpr = exprs (data[[1]]) will just give an object with one header line (      GSM248431 GSM248432 GSM248433 GSM248434 GSM248435 GSM248436 GSM248437 GSM248438 GSM248439 GSM248440 GSM248441 GSM248442)

Well, if you use a SOFT formatted files:
#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that "Raw data provided as supplementary file. Processed data is available on Series record". So one may want to check series record instead:
# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz)

# and check what´s inside:
exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?
If you want to download them, you can:
filePaths = getGEOSuppFiles("GSE9838")
filePaths
untar("GSE9838/GSE9838_RAW.tar", exdir = "GSE9838/CEL")
list.files("GSE9838/CEL")

...what will allow you to see all the CEL files with the data.

Thank you so much! I will try to open the raw data.


text:
scaling algorithm of deeptools
I am using deeptools to generate a metagene plot. I am wondering how the deeptools scale genes if the size of genes is shorter than bin size. For example, if my bin size is 1000 bp and the gene is about 450bp, how does the deeptools computeMatrix scale it to make a metagene plot of 1000 bp?


text:
filter secondary alignments using pysam
I am new to python and trying to learn. The below is an attempt to filter out secondary reads in a bam using pysam. I ran mark duplicates and got an error on several bam files due to secondary alignments. I added comments to each line to illustrate my thought process.  Thank you :).
The specific mark duplicates error:
Value was put into PairInfoMap more than once

Here is the script:
#! /usr/bin/python  ## call python script
import sys ## import python system functions
import pysam  ## import module

bam = pysam.AlignmentFile(*.bam, "rb")  ## open bam and read
click.echo("Reading BAM file")  ## output message
for read in bam:  ## start loop and iterate over each bam
    if read.is_secondary=true  ## not the primary alignment
    if read.has_tag('XS')  ## look for Secondary Alignment Score tag
        XS = read.get_tag('XS')  ## read and store AS value
    if read.has_tag('AS')  ## look for Alignment Score tag
        AS = read.get_tag('AS')  ## read and store value os XS tag
    if AS < XS  ## Alignment score less then
        read.write(read)  ## only print primary alignments
bam.close()  ## end processing

Is there an issue with the script?  If so what is it?  Are you looking for feedback on the script?

Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') continue block.
I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.
for read in bam.fetch()
    if read.is_secondary=true
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')
        continue
    if read.has_tag('AS')
        AS = read.get_tag('AS')
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)
bam.close()

I have not tried to executec it as this is my first python script i am looking for suggestions. Thank you :).

I would suggest running the script to determine if it is doing what you want.  What I pointed out above is def a flaw in the logic.

I ran the script (modified it in the post) and there are many syntax errors... being new to python I am not sure what is best... what do you suggest? Thanks :).

I would start will researching how to write to a file.  read.write(read) is not correct syntax.  You will need to open a file and write to that.  Go error by error until you have a working script then determine if the out is correct.


text:
CDS length for each human gene
Does anyone know where and how could I download a list of all human genes and the length of the coding sequence for each gene? Is it possible to do this on the NCBI site, ensembl?

Which coding sequence? I mean, do you just want whichever has been designated the 'canonical' transcript or do you want all possible isoforms?

Hi terdon, thanks for the quick reply! Yes exactly, the canonical transcript is good enough!

Michael G. suggests to take a look at relevant front-end, NCBI's eFetch. Which is supposedly perfect for what you need.

While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them).
Essentially, you just need to go to BioMart, and

select "Ensembl Genes 96" (the number will change if the version changes) as the database and "uman Genes" as the dataset.
Click on "Filters", and set Gene type to coding and Transcript type to protein_coding.
From "Attributes", select whatever you want to see. The "CDS Length" is under "Structures".

Go to the NCBI homepage and search for 'human genome' (the dropdown menu on the left-hand side of the search box should be 'All Databases'). In the result page, click on the 'Download' button; choose 'RefSeq' as the 'Source database' and 'Feature table' as the 'File type' as shown in the screenshot below:

This will download a genome_assemblies.tar file which contains the feature table file, GCF_000001405.39_GRCh38.p13_feature_table.txt.gz. This file has all of the annotated features in a tab-delimited format with the following fields:
              # feature [  1]: gene
                  class [  2]: transcribed_pseudogene
               assembly [  3]: GCF_000001405.39
          assembly_unit [  4]: Primary Assembly
               seq_type [  5]: chromosome
             chromosome [  6]: 1
      genomic_accession [  7]: NC_000001.11
                  start [  8]: 11874
                    end [  9]: 14409
                 strand [ 10]: +
      product_accession [ 11]:
   non-redundant_refseq [ 12]:
      related_accession [ 13]:
                   name [ 14]: DEAD/H-box helicase 11 like 1
                 symbol [ 15]: DDX11L1
                 GeneID [ 16]: 100287102
              locus_tag [ 17]:
feature_interval_length [ 18]: 2536
         product_length [ 19]:
             attributes [ 20]: pseudo

You can parse this file to obtain the CDS lengths of all annotated proteins. As an example, you can get a list of all annotated RefSeq proteins and their lengths as follows:
zcat GCF_000001405.39_GRCh38.p13_feature_table.txt.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($1~/^#/ || $1=="CDS")' \
    | cut -f6-11,13,15,16,18,19

The output of the above command will have the following fields:
             chromosome [  1]: 1
      genomic_accession [  2]: NC_000001.11
                  start [  3]: 69091
                    end [  4]: 70008
                 strand [  5]: +
      product_accession [  6]: NP_001005484.1
      related_accession [  7]: NM_001005484.1
                 symbol [  8]: OR4F5
                 GeneID [  9]: 79501
feature_interval_length [ 10]: 918
         product_length [ 11]: 305

Ensembl has an FTP site that allows you to select and download only the coding sequences from many different genomes.
https://useast.ensembl.org/info/data/ftp/index.html
To determine the length of those sequences, download the associated gtf or gff3 annotation file. The annotation file is tab delim. The fourth and fifth column represent genomic loci of the annotated region. Subtract the amount in the fourth column from the amount in the fifth column to yield the length of all the annotated features.
You can easily do this in R after loading the file in the environment using the Magrittr library. The following code will create a new column called gene_length with the associated gene lengths.
install.packages("magrittr")
# this only needs to be done once
library(magrittr)
# must be run each time the library is neaded
annotation.gtf <- read.table("path/to/annotation.gtf")
annotation.gtf$start <- annotation.gtf[,4]
annotation.gtf$end <- annotation.gtf[,5]
annotation.new-column.gtf <- annotation.gtf %>%
          mutate(gene_length=end-start)

authors: Александр Глебус


text:
How can I calculate coverage at single bases using a bam file?
I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:
chr base_position total_coverage fwd_coverage rev_coverage

I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information.

Note: not yet tested, so there may be some additional fiddling with command line options needed
The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions):
samtools depth -a in1.bam > depth_in1_both.tsv

To split this by forward and reverse, you can use an initial pipe through samtools view to exclude or include reverse-complement mappings:
samtools view -b -F 0x10 in1.bam | samtools depth -a > depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a > depth_in1_rev.tsv

To filter for particular regions, this can be pre-piped through bedtools intersect:
bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a > depth_in1_regions_fwd.tsv

Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...
If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using samtools depth -aa), then the lines can be pasted, with anonymous pipes, with unnecessary fields filtered out using cut:
paste <(samtools depth -aa in1.bam) \
      <(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      <(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

Or equivalently, if using intermediate files that have been generated already:
paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

It would also be possible to use join to combine the files (in cases where file output is not consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.

git clone https://github.com/lh3/htsbox
cd htsbox && make
./htsbox pileup -cCf ref.fa aln.bam | less -S

This output a VCF containing positions covered by at least one read. Tag ADF and ADR give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.

Brent Pedersen (@brentp) claims mosdepth is 2x as fast as samtools depth. I'm not sure whether this provides strand information, but it might be the fastest tool available.
Another tool of possible interest (also written by Brent) is indexcov. This doesn't give single bp precision, but is able to provide even quicker estimates of coverage at 16Kb resolution.
Installing these tools from source is non-trivial, and I struggled with installation in the past. However, they are both now available in bioconda, which should make installation simple.


text:
Error in running Mutect2
I am trying to run Mutect2 on a .bam file but I get an error by Googling can not be tackled
Have you seen this error?
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ ./gatk Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
Using GATK jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar Mutect2 -R reference.fa -I sample.bam --germline-resource af-only-gnomad.hg38.vcf.gz --panel-of-normals pon.vcf.gz -O single_sample.vcf.gz
13:16:43.616 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/data2/RNASeq/Angel/gatk/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
Jan 31, 2022 1:16:43 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
13:16:43.733 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - The Genome Analysis Toolkit (GATK) v4.2.4.1
13:16:43.734 INFO  Mutect2 - For support and documentation go to https://software.broadinstitute.org/gatk/
13:16:43.734 INFO  Mutect2 - Executing as fi1d18@RBGO-Server2 on Linux v4.15.0-166-generic amd64
13:16:43.734 INFO  Mutect2 - Java runtime: OpenJDK 64-Bit Server VM v11.0.13+8-Ubuntu-0ubuntu1.18.04
13:16:43.734 INFO  Mutect2 - Start Date/Time: 31 January 2022 at 13:16:43 GMT
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.734 INFO  Mutect2 - ------------------------------------------------------------
13:16:43.735 INFO  Mutect2 - HTSJDK Version: 2.24.1
13:16:43.735 INFO  Mutect2 - Picard Version: 2.25.4
13:16:43.735 INFO  Mutect2 - Built for Spark Version: 2.4.5
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.COMPRESSION_LEVEL : 2
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
13:16:43.735 INFO  Mutect2 - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
13:16:43.735 INFO  Mutect2 - Deflater: IntelDeflater
13:16:43.735 INFO  Mutect2 - Inflater: IntelInflater
13:16:43.735 INFO  Mutect2 - GCS max retries/reopens: 20
13:16:43.735 INFO  Mutect2 - Requester pays: disabled
13:16:43.735 INFO  Mutect2 - Initializing engine
13:16:44.015 INFO  Mutect2 - Shutting down engine
[31 January 2022 at 13:16:44 GMT] org.broadinstitute.hellbender.tools.walkers.mutect.Mutect2 done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2155872256
***********************************************************************

A USER ERROR has occurred: Cannot read file:///data2/RNASeq/Angel/gatk/pon.vcf.gz because no suitable codecs found

***********************************************************************
Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace.

I also get this
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/gatk$ file pon.vcf.gz
pon.vcf.gz: HTML document, ASCII text, with very long lines

Sorry, I think the nature of stackexchange is such that you will need to post a new question to get help with the new problem. Perhaps a better location for this question is biostars? https://www.biostars.org/

Yep, edits are for clarification, questions can't be moving targets (I understand it's frustrating, but changing the question renders the answers obsolete and that just generates a mess - which is the reason we don't allow it here)

Your input file pon.vcf.gz should be a gzip file, but from your file command it looks like it is an HTML file instead. Perhaps you downloaded the VCF file incorrectly?
Here is what I get on a bona fide gzip file in Ubuntu:
$ echo hi | gzip -c > hi.gz
$ file hi.gz
hi.gz: gzip compressed data, from Unix, original size modulo 2^32 3

Thank you, you were right, I downloaded the proper file but I get the same error, please have a look at my main post]

Downloading that vcf file from this link fixes the problem (but generates a new one)
wget https://storage.googleapis.com/gatk-best-practices/somatic-hg38/1000g_pon.hg38.vcf.gz


text:
An error with annovar
I want to install annovar but in one particular part I always get this error
i1d18@RBGO-Server2:/data2/RNASeq/Angel/annovar$ perl annotate_variation.pl -downdb -webfrom annovar avsnp147 humandb/
NOTICE: The --buildver is set as 'hg18' by default
NOTICE: Web-based checking to see whether ANNOVAR new version is available ... Done
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.gz ... Failed
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.idx.gz ... Failed
WARNING: Some files cannot be downloaded, including http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.gz, http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.idx.gz
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/annovar$

Googling did not return any clue
Any help please?

My first attempt would be to change the version of the genome build to something more recent than hg18.

I don't think hg18 is the issue, according to their documentation hg18 is supported. However, the commands you used are not exactly the same as they showed in the tutorial. Did you try to do exactly what did suggested? If that still does not work, I would suggest to open an issue on their github.

Changing to hg38 solved the problem, thanks


text:
Proper use of BWA MEM on multiplexed GBS sample
I have a multiplexed lane of GBS sequencing reads as a fastq file. I understand the first step is to demultiplex and trim the adapter sequences from the reads.
This yields many individual fastq files that correspond to individual libraries.
I would like to genotype my samples, but I need to map the reads first. My question is how should I perform the BWA MEM alignment?
Do I perform BWA MEM on each of the resulting fastq files that correspond to individual libraries...
OR
I should I perform BWA MEM on the large original fastq file that contains all my sequence reads before demux and read trimming?
Forgive my ignorance, I've only ever used TASSEL (which works great) and I'm trying to expand my understanding and flexibility.

The first thing to do is to group the FastQ files together following this post. Once you have your reads grouped, you could trim your files before mapping. Then you could perform your BWA alignment. Since you mentioned MEM, I suppose your read length should be between 70-100bp as suggested by the author:
It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for
longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such
as long-read support and split alignment, but BWA-MEM, which is the latest,
is generally recommended for high-quality queries as it is faster and more accurate.
BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

Then you could run your BWA MEM, supposing you have 10 CPUs available, like so:
bwa mem \                                  # bwa mem algorithm
-M \                                       # mark shorter split hits as secondary (for Picard compatibility).
-t 10 \                                    # CPUs thread
reference.fasta \                          # reference genome to be aligned to
< sample_r1.fastq \                        # fastQ read 1 (can be .gz)
< sample_r2.fastq | \                      # optional FastQ read 2 (can be .gz)
samtools view -bS -@ 10 > sample.bam       # pipe to samtools to convert SAM to BAM on the fly.

Hope this helps !

Trimming is something you always want to do before mapping to make sure that the bases at each position are of sufficient quality, and to make sure that you are not aligning your adapters to the reference genome.
I am not entirely sure what the biological difference is represented in your individual libraries (that depends on your research question). But most likely you want to demultiplex the reads before mapping with BWA as well. As a result you would indeed get a sam/bam file per library. Let me know if anything is still unclear!


text:
Need help in identifying origin of variants in CIVIC
We usually denote the origin of a mutation as either somatic or germline. This information is usually available in certain databases such as CIVIC, ClinVar, COSMIC etc. But when we come to variants that have been denoted as "Overexpression", "Upregulation" , "Methylation" etc., what would be the origin? Any leads would be very helpful. Thanks in advance.
For example from Cristina Aguirre‐Portolés et al., 2018:

ABCA1 is significantly overexpressed in patients at advanced stages of
  colorectal cancer, and its overexpression confers proliferative
  advantages together with caveolin‐1 dependent‐increased migratory and
  invasive capacities. Overexpression of ABCA1 was associated with
  overall survival independently of the molecular subtype of the cancer.

Does this overexpression phenomenon arise in the germline or does it happen somatically? How do we conclude?

I don't really think you can if it wasn't mentioned in the paper...

@terdon Does it make sense that when we say "Overexpression", we refer to expression of a gene using normal expression as control? As in, when we measure the expression level we also take a normal from healthy tissue. So doesn't it mean that this overexpression is in fact a somatic event?

There is no one answer to that. It will always depend on the specific context. The control could be normal tissue vs tumor, yes, but it could also be healthy individuals vs patients. Or 10 days after diagnosis vs. 2 years after diagnosis. In the example you cite, we cannot know if the overexpression of ABCA1 is limited to the tumor cells or systemic in the entire body. Perhaps there is an immune response to the tumor that stimulates systemic overexpression. There could be any number of possibilities.

Actually, it does seem likely that they are referring to the tumor since they mention increased migratory and invasive capacities and Overexpression [...] was associated with overall survival independently of the molecular subtype of the cancer, but it isn't clear from your extract what they are comparing to. In any case, there will not be one single answer here and it will always depend on the specific study.

@terdon Ok that makes it a lot clearer now. Thank you for your help !


text:
ANNOVAR, 1000g database calling
I have recently downloaded Annovar, and succesfully filtered my files with a few databases, downloaded from Annovar. But when I try to do the exact same thing with the 1000 genomes project database, I get this error many times in the terminal, one for each line of **.
Argument "&" isn't numeric in numeric eq (==) at annotate_variation.pl line 1962, <DB> line &&.

& can be A, C, G, or T.
&& every line of the database
I have tried filtering with 1000g, 1000g2010 and 1000g2010jul databases, but all give these errors.
An example of the commands I use for filtering:
perl annotate_variation.pl --filter --dbtype CEU.sites.2010_07 testfile1.txt humandb/

Can anyone tell me if I'm doing something wrong, if this is known, and if I can do anything about it?

Is it an error (i.e. it stops the script from proceeding further) or a warning (i.e. alerts you with a message but doesn't stop execution)? IIRC this is a common warning with most ANNOVAR operations.

That's a warning. Presumably, a match failed ($& is the default Perl variable which holds the result of a match operation). So something in your input is not exactly what the script was expecting. You can probably safely ignore these messages.
Alternatively, update your question and add a few example lines of input which we can use to reproduce the error and maybe we can figure something out. If that isn't possible, at least show us line 1962 of the annotate_variation.pl script since that's the one that's failing.


text:
What are values of FILTER column of vcf files produced by Mutect2
I have called SNVs with Mutect2, now in filter column of a file called filtered vcf I have a lot of things like
> length(unique(filtered.vcf$FILTER)) [1] 409 >

> unique(filtered.vcf $FILTER)
  [1] "germline"
  [2] "PASS"
  [3] "clustered_events"
  [4] "clustered_events;panel_of_normals"
  [5] "weak_evidence"
  [6] "haplotype"
  [7] "base_qual;haplotype"
  [8] "clustered_events;germline;haplotype"
  [9] "base_qual;clustered_events;germline;haplotype"
 [10] "clustered_events;haplotype"
 [11] "clustered_events;haplotype;weak_evidence"
 [12] "base_qual;clustered_events;haplotype"
 [13] "haplotype;map_qual"
 [14] "germline;panel_of_normals"
 [15] "base_qual"
 [16] "germline;haplotype"
 [17] "germline;panel_of_normals;position"
 [18] "clustered_events;orientation"
 [19] "clustered_events;orientation;strand_bias"
 [20] "base_qual;haplotype;map_qual"
 [21] "clustered_events;fragment;haplotype"
 [22] "base_qual;clustered_events"
 [23] "germline;haplotype;panel_of_normals"
 [24] "clustered_events;germline;haplotype;panel_of_normals"
 [25] "fragment"
 [26] "clustered_events;haplotype;map_qual"
 [27] "clustered_events;germline"
 [28] "clustered_events;germline;panel_of_normals"
 [29] "germline;weak_evidence"
 [30] "orientation"

I tought that the filtered.vcf Mutect2 output file contain only the filtered events, but the values I see suggest otherwise. How should I understand the output?

You haven't told us what commands you ran, so I am assuming you first ran Mutect2 and then FilterMutectCalls. If so, the FilterMutectCalls does not remove any variants from its input, it will simply flag any that fail its builtin filters. Any variants that did not fail a filter are marked as PASS and the rest have a filter added.
Now, whether you want to filter these out or not is up to you. You need to look at the VCF header where each of the filters is defined (admittedly, not very clearly) and make a judgment call on whether you consider this filter to be a deal breaker or not.
For instance, the germline flag means:
##FILTER=<ID=germline,Description="Evidence indicates this site is germline, not somatic">

So that indicates that gatk believes this is more likely to be a germline variant than a somatic one, that it is more likely to be a constituent variant present in the sample's germline cells than a variant that is only found in the sample's tumor cells.
Another clear case is the base_qual filter which means:
##FILTER=<ID=base_qual,Description="alt median base quality">

When a variant is flagged with this filter, it means that the median fastq base quality of the reads supporting the variant was significantly lower than the median fastq base quality of the reads supporting the reference alele at the same position, making it likely that the variant is a false positive.
You need to read the VCF headers and any gatk documentation you can find (warning: these filters are not very well documented at all, in my experience), understand what the filters are and then decide what variants you consider real based on what you know about your sample, your experimental design and the question you are trying to answer.

How did you get to the explanation of base_qual from just the "alt median base quality"? I'd like to explain each filter in simple English like you've done for base_qual above and there seems to be very little proper documentation.

@RamRS yes, extremely little. I got this one from parsing  "alt median base quality" as "median base quality in the reads supporting the alternative allele" and it has also been borne out by manual testing. Unfortunately, there is almost no documentation available at all for these things.


text:
How to make a violin plot around quasirandom dots
I want to have a plot violin plot around quasirandom dots - like below

This is my data
> str(long1)
Classes ‘data.table’ and 'data.frame':  116 obs. of  4 variables:
 $ gene    : chr  "WT1" "KRT7" "MUC16" "EPCAM" ...
 $ variable: chr  "1016" "1016" "1016" "1016" ...
 $ value   : num  5.82 8.41 11.22 7.21 5.4 ...
 $ group   : chr  "Tumour" "Tumour" "Tumour" "Tumour" ...
 - attr(*, ".internal.selfref")=<externalptr>
>

By this code I get plot which does not have any violin
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

I also tried adding geom_violin like
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_violin()+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

But gives a weird plot not similar to what I want

Any help please
EDITED
Thanks to @haci
By this code I get this plot but as you see the point color inside the violin is not really similar to what I want
Any suggestion to get a similar color scheme for dots inside the violins
p <- ggplot(long1, aes(x=group, y=value,fill=group))+guides(fill = "none")
p +geom_violin(alpha = 0.1)+scale_fill_manual(values=c("#fcf4e0", "#fcf4e0"))+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,size=2,stroke=0.5,aes(colour = variable))+scale_color_hue(l = 5, c = 150)+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right",legend.text=element_text(size=12)) +
  theme(legend.title=element_blank(),legend.text.align = 0,text = element_text(family = "Arial"),
        axis.title.x = element_text(size = 0.00001),
        axis.title.y = element_text(size = 17))+theme(axis.text = element_text(
          color="black",
          size=12,
          face=1)
        )+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 17,face="italic"))+ theme(axis.line = element_line(colour = 'black', size = 0.8))

You are almost there, if you set the alpha parameter within your geom_violin() call (to sth like 0.1) you should be able to get a plot similar to the example plot.
The fact that you only have two data points for the Stroma group and these two points having similar values, you will end up getting a horizontal-line like violin for this group. Moreover, you do not have that many data points in general anyway so you might consider not using violins at all, especially if the outliers in Tumour group is coming from the same (maybe somewhat problematic?) sample, that might need to be dropped from the analysis. The reader should be able to grasp the essence of the data with our without violins (violins could help making the outliers obvious though).

Thanks @haci . In edited version of my post I am struggling with the same color scheme, please do you have any suggestion here?

What you mean @M_

@Exhausted Don't edit your question and branch off, ask a new question as a follow up (by mentioning this question) instead. When you edit questions that have been answered, you confuse future users by removing context. Scientific forums are not your personal help desk.

@Exhausted, you can try increasing the alpha parameter in your geom_quasirandom() call (or remove that parameter altogether).


text:
A strange error while clustering of single cell RNA-seq using URD package
I am using URD R package on my single cell RNA-seq data hoping to be able to trace gene expression changes during the development. After dimension reduction step by PCA and diffusion map whatever I am trying to cluster some cells, but whatever I do I am getting this error
> object <- graphClustering(object, num.nn = c(5,8,10,15), method="Louvain", do.jaccard = TRUE)
Error in graphClustering(object, num.nn = c(5, 8, 10, 15), method = "Louvain",  :
  unused arguments (num.nn = c(5, 8, 10, 15), method = "Louvain", do.jaccard = TRUE)

> print(object.6s)
URD object: 12451 genes x 1556 cells.

object.6s.mnn <- graphClustering(object.6s)
Show Traceback Rerun with Debug Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘graphClustering’ for signature ‘"URD"’

I tried to change Louvain, jaccard, ... but I am still getting this error. How can I get rid of it?

The error message says that you should remove the num.nn, method and do.jaccard argument. Did you try it?

Thank you, yes I removed all except for my URD object but the error is permanent :(

Can you post the new command and the errors? Could you add the 'print(object)' result?

Please [edit] the question to add the code (in code formatting) . It seems like the error is different now. I'll research a bit

The problem is they define a new class but they don't define methods for said class. The package would need (IMHO) more test and add methods for said class. It would be still possible to work with functions if they defined a generic class, but that's beyond what we can answer/solve here, they need to modify the package.

Also, please consider asking to the right people next time. If you post an issue in the package it signals that there is something wrong with the package(you followed the tutorial as you did and you found a problem), while if you post it here I understand that there is some problem while using the package (you are trying to do something not covered by the help file or the tutorial and don't know how to do it). Either way if it is a "simple" informatics error probably Bioinformatics.SE is not the best place to post it.

Thank you I also posted in GitHub but URD package developers did to reply

Yes, I saw, it but don't expect people to answer the emails or questions in 5 minutes. The developers (and people on the forums, here on biostars or in support.bioconductor) have other things to do and are answering you on their own time, because they want to help you. To help you they need to spend some time figuring the error and correcting it for you, so it takes some time to give you an answer. I hope you get your answers and solution to these problem.

definitely you are alright, I am also thankful for your time

This solved the error
object <- URD::graphClustering(object)

By this function I am telling R to use the graphClustering() function from the URD namespace. I probably had loaded another package with a different graphClustering() function - or defined one myself - which is taking precedence over URD's function.

It's nice that you have shared a solution to your problem. But since it was not very clear what is the actual problem, this answer is not really explaining anything.

I added the possible reason


text:
predict the IC50 after using Schrodinger suite Glide tool for docking
How to predict the half-maximal inhibitory concentration (IC50) of my ligands.
I'm currently using Schrodinger suite software and I don't know how to predict the IC50?

It is usually helpful to do a little research with a search engine or to give more information about what exactly you need. See e.g. https://www.schrodinger.com/kb/936

I don't think it is possible just from the docking scores. Docking scores are designed to give you an estimate of the binding affinity (entalphy), but they do not take entropic contributions into account very well. The best you can do are molecular dynamics simulations to calculate the free energy of binding. Which can be computationally quite costly e.g. with GROMACS or CHARMM.
You can estimate the IC50 with molecular dynamics simulations as well, but that is not straightforward. E.g. with Umbrella simulations. An example for a potassium channel and a very simple ligand can be found here.
Another approach would be to build a predictive model. If you have a library of compounds at hand you could use the chemical structures and e.g. the docking scores to estimate the IC50 based on statistics or machine learning. But biology is very messy in general, though if you get a representative dataset with very similar compounds to the one that you have in your dataset then you might be successful. You have to be sure that the binding site and pose are correct though.


text:
How To display bar graph with error bars?
I am using the below code to get the bar plot. I want to get the bar plot with error bars displayed.
import seaborn as sns

def plot_coefficients(target, coeffs_df, n=10):
    if coeffs_df.shape[0]==0:
        return None # nothing to plot here!
    data = coeffs_df.copy()
    data = data.sort_values(('Coefficient','mean'), ascending=False).head(n)
    g = sns.barplot(
        data=data,
        y='Gene',
        x=('Coefficient','mean'),
        hue='Feature Type',
        hue_order=[
            h for h in [
                'Genomic Variant',
                'Gene Expression',
                'Proteomics',
                'Methylation',
            ]
        ],
        dodge=False
    )
    g.set(
        xlim=(0,1),
        xlabel='Importance [0,1]',
        ylabel='',
        title=f"Top {target} predictor features"
    )
    return g

Coefficient here is a score generated from predictive modelling that ranges from 0 to 1. 2. data here is entire data and I created a mean column. Ideally I wanted to calculate it in seaborn plot rather than creating a mean column separately before.

I presume you import seaborn as sns, right?

@KamilSJaron Yes

I did put forward a comment, it was removed I think the idea is its a full post or zero. The problem proving a negative is tricky.
I don't think its doable within the exact system specified, i.e. a Seaborn bar chart. Again I could be wrong.
If the OP went to a Seaborn dot line formation - that would work sns.barplot is changed to sns.pointplot but then all the nice colouring disappears. Standard deviation  is one error range if all the values are present: but I think the OP will need to request it (.sd() ?? can't remember).
In matplotlib there is an error bar method where the coder specifies the error bars and that is plt.errorbar. An example from the official docs is here. Thus call a bar chart and then add the error bars onto it.
Summary The reasons I don't think its doable in Seaborn:

I'm pretty certain Seaborn has no independent errorbar method unlike matplotlib
Seaborn can do error bars given the full data set, but those are created automatically for a given method. Maybe there is a "sns.barplot-errorbar" method but I don't know it.

The Seaborn approach does carry the advantage of not having to work out the error bar independently.


text:
Import gene list to Seurat to define cell types
Is there a way to import gene list into Seurat to define cell type? The default cell types in Seurat is not enough for our research. For example, we want to mark a subtype of B cells in Seurat, but seems like Seurat only have general B cell, I'm planning to download gene list from Genomic Cytometry, and then import the gene list into Seurat and then Seurat define cell types in the cluster based on the gene list I import. Does Seurat support this function?

Welcome to the site. Did you find anything on the vignette and on the several questions (and answers) posted online?

I checked, but unfortunately I couldn't find any information.

I am in roughly the same situation: I have a list of marker genes for each cell type, and I now I want to assign a cell type to each cell using this. Did you ended up finding solution and would you mind sharing?

You can use scID with your gene list to find the cells that match your gene list. You can then find which cluster in Seurat have these cells.

Seurat does not define cell types by name. It clusters and assigns each cell to a cluster, from 0 to X.
If your data has the cell type (e.g. B,T, Mast cells) it means that someone annotate the clusters so that they have a biological meaning.
You can assign different names to the clusters by using the AddMetaData function.
sampleID  cellType
cell1     B
cell2     B
cell3     T
cell4     T

newMeta <- oldMEta
newMeta$origType <- newMeta$cellType
newMeta$cellType <- c("B1","B","Treg","Tkill")

object <- AddMetaData(object, newMeta)

The new metadata should look like
sampleID  origType  cellType
cell1     B         B1
cell2     B         B
cell3     T         Treg
cell4     T         Tkill

You can then swap the labels using SetAllIdent(object = test, id = "cellType"). Your tSNE or uMAp, will now show the new cell types.
Have a look at Manually define clusters in Seurat and determine marker genes

Thank you so much for your answer. But actually I want to define a new subtype of B cells in the B cell cluster. I know the genes that express the new subtype of B cells, but the current Seurat cluster only showing the general B cell. I'm thinking about doing a subtype of B cells, and then import the the gene list to the Seurat to see if Seurat can show the subtype of B cells based on the gene list I provided. I'm not sure if Seurat can do this.

Seurat does not use gene lists to cluster the data. What I usually do to look for sub-cell types, is to subset the dataset (only the cells from one cluster) and re-run the clustering/PCA steps to see how those cells (and only those cells) are clustering together. The expression of your genes will tell you if the sub-clusters contains specific sub-types.  Have a look at the Seurat documentation under "Selecting particular cells and subsetting the Seurat object".

Hello, you mentioned that Seurat not use gene list to cluster data, but I thought that Seurat cluster cells based on the expression of certain genes, if a cells have high expression in certain genes, they will be define as a specific cell types, do you mean that this is not the case?

Seurat clusters based on the expression of certain genes, but those genes will be determined independently for each experiment. By default, it would be the most variable genes. It does not know if those genes have been previously described to be associated with a specific cell type.

Maybe you can try Seurat::AddModuleScore(), then FeaturePlot() and see if some of your B cells are different. After plotting this on GenePlot(), perhaps you can set a cutoff, then assign identities.
Alternatively, use your B cell gene list in RunPCA(object, pc.genes = yourgenelist) instead of the usual variable genes. You can also subset B cells first with SubsetData(), and run the above.