What does an output of -I for all amino acids mean in a psi-blast pssm matrix file?

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


The prior count ends up getting scaled by the ratio of a library size to the average library size and then multiplied by 2 before getting added to each library size (I'm sure there's a good reason for that, but I don't know what it is). Using the example df data frame in your post, let's walk step-by-step through what cpm() is doing:

df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)
prior.count = 0.5
# First, we need to calculate a library size
lib.size = colSums(df)
# Calculate the average library size and the adjusted priors
ave.lib = mean(lib.size)
adjusted.prior = prior.count * lib.size / ave.lib
# Update the library sizes
adjusted.lib.size = lib.size + 2*adjusted.prior
# Now we can compute the CPM
customCPM = t((log(t(df) + adjusted.prior) - log(adjusted.lib.size) + log(1000000))/log(2))

The matrix transposition looks really ugly in the last step, but it's needed so R adds things like adjusted.prior to each row of df.

Note that the prior count is only used if you compute log(CPM), it's unused otherwise.


For sharing, you can test the proportion of overlapping peaks using bedtools intersect or the find.overlap function in GenomicRanges packages in R.
What makes more sense is perhaps to merge all the peaks together using bedtools merge, and then count the number of reads/fragments in each peak for each sample using featureCounts. The output of this is a matrix, which is very similar to gene expression results. Then you can do some differential tests on this data.


How can I get bedtools to tell me which genes are being expressed?

Bedtools is not the best tool for this. You can use featureCounts, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.
bedtools coverage is used to calculate the overlap between each read in BAM and each line in BED. It is not very straightforward to quantify gene expression. Plus bedtools does not account for exon in each gene and summarize them to gene level, like featureCounts does.


How to deal with duplicate genes having different expression values?

I assume you're familiar with the various issues surrounding FPKMs, so I'll not expound upon them.

As a general rule, you should be using gene IDs rather than gene names, since the former are unique while the latter are not. If you only have access to data quantified on gene names, then the appropriate way to merge RPKMs is with a weighted sum:

$FPKM_{gene} = \frac{FPKM_{copy1} * Length_{copy1} + FPKM_{copy2} * Length_{copy2}}{Length_{copy1} + Length_{copy2}}$

As an aside, rRNA expression levels will tend to be wrong, since they're normally excluded by either poly-A enrichment or the use of ribo-zero.


What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?

First, I'll observe that Freeman centrality as you define it is the reciprocal of the Barycenter centrality as you define. It looks like people often use the reciprocal instead of the original definition, so in that sense they do look the same.
I found the following comparison here (though it doesn't name it as Freeman):

There are 2 types of distance centrality scores, Closeness Centrality
and Barycenter Centrality. Barycenter Centrality for vertex v defined
as:
1 / (total distance from v to all other vertices)
Closeness scores are calculated using the formula 1 / (average
distance from vertex v to all other vertices) and Barycenter scores
are calculated as 1 / (total distance from vertex v to all other
vertices).

So it appears that for this package it is a difference of the reciprocal of the total vs. reciprocal of the average.
Notably, wikipedia differs here, in that it matches your source:

Freeman's closeness centrality, the total geodesic distance from a given vertex to all other vertices, is the best known example.citation

Yet another source notes the existence of a "normalized" version of Freeman's closeness centrality, which normalizes centrality of a node by the degree of the network. It also comments on whether it is the reciprocal or not:

Closeness (Freeman '79)
The graph-theoretic distance of a given node to all other nodes. The sum of the rows/columns of the geodesic distance matrix of a graph.
Simple closeness is an inverse measure of centrality: the larger the numbers, the more distant an actor is, and the less central. Should really be called "farness".
Normalized version divides the minimum "farness" possible (N-1) by "farness" to simultaneously make the range 0 to 1 and invert the measure so that larger values correspond to greater centrality (truly "closeness")

So I guess that it might depend on what tools you are using. Most of the time, it might be the same. Other times, it might be normalized in different ways with either averages or for 0-1 ranges.
Closeness centrality is probably the more standard term, so it is unlikely that anyone would criticize it. But it appears that there is sometimes some variation in what people mean by that term and whether they are using a normalized or a raw closeness centrality.
In my opinion the normalization leading to a range from 0-1 is nice because it makes the measure at least somewhat generalizable across networks. But it's up to what you need to do.


Are there genomics databases for neurological and psychiatric disorders?

Try with PsyGeNet:

  PsyGeNET (Psychiatric disorders Gene association NETwork) is a resource for the exploratory analysis of psychiatric diseases and their associated genes

It is based on literature mining and curated by experts.

Here is the publication if you want more details:

PsyGeNET: a knowledge platform on psychiatric disorders and their genes

and they even have a Bioconductor package for programatic queries:

psygenet2r: a R/Bioconductor package for the analysis of psychiatric disease genes


Is there a straightforward way to get mismatches/indels in a BAM file using pysam?

You're looking for pileup, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            ...handle mismatches...


You could give a try with dplyr's full_join() and purrr's reduce(). full_join() will make sure that you won't lose any data and reduce() will merge tables in your list from left to right, left being the first table of the list in the first iteration but would become merged table1+table2 in the second iteration and so on. This left to right "reduction" will also be helpful when renaming the columns of the resulting table.

Without your workspace I have attempted to replicate your data: There is a list of 3 data frames, each with 2 columns:

library(dplyr)
library(purrr)

df1 <- data.frame(CDR3.aa = c("A", "B", "C"),
                  Count = c(10,20,100))

df2 <- data.frame(CDR3.aa = c("A", "B", "D"),
                  Count = c(100,200,300))

df3 <- data.frame(CDR3.aa = c("D", "E", "F"),
                  Count = c(10,20,100))

my_list <- list(df1 = df1, df2 = df2, df3 = df3)

> my_list
$df1
  CDR3.aa Count
1       A    10
2       B    20
3       C   100

$df2
  CDR3.aa Count
1       A   100
2       B   200
3       D   300

$df3
  CDR3.aa Count
1       D    10
2       E    20
3       F   100

# "by" argument specifies which column to use when merging
my_table <- reduce(my_list, full_join, by = "CDR3.aa")

names(my_table) <- c("CDR3.aa", names(my_list))

> my_table
  CDR3.aa df1 df2 df3
1       A  10 100  NA
2       B  20 200  NA
3       C 100  NA  NA
4       D  NA 300  10
5       E  NA  NA  20
6       F  NA  NA 100

Clonotypes that are not present in table will be "NA".


I personally find the GDC portal not very user-friendly. However, some of the data is available from other sources. For example, you can download TCGA-BRCA FPKMs and phenotypes from UCSC Xena here: https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Breast%20Cancer%20(BRCA)

All the sample info is already combined and the labels are the classic TCGA barcodes.


How to get the corrected matrix after SVA batch effect correction

If you want to plot the "corrected" expression, you will need to remove the variation introduced by these surrogate variables. Removing the expression affected can introduce some bias too and it is usually not recommended (despite comBat doing so). You should apply linear algebra, you can look at here is an example how to do it:

library("corpcor")
s <- corpcor::fast.svd(t(scale(t(dat), center = TRUE, scale = TRUE)))
pcSds <- s$d
pcSds[1] <- 0
svdexp <- s$u %*% diag(pcSds) %*% t(s$v)
colnames(svdexp) <- colnames(dat)

But if your surrogate variables do have a big enough influence on your MDS plot to be observable (IMHO) there would be too much batch effect to work with the data.


What you are looking for is facetting. Using this key word on any search engines you will find dozen of answers describing its use: https://ggplot2.tidyverse.org/reference/facet_grid.html and https://ggplot2.tidyverse.org/reference/facet_wrap.html

To prepare the faceting in ggplot2, a possible way is to bind both dataframes together by first specifying which group they belong such as:

cytokine$Group = "Cytokine"
il17$Group = "IL17"
DF <- rbind(cytokine,il17)
DF

    gene       variable value    Group
1:  CCT5     Res_immune 10.60 Cytokine
2:  CCT5 Non_res_immune 10.50 Cytokine
3:  CDK1     Res_immune  9.25 Cytokine
4: CXCL5 Non_res_immune  9.82 Cytokine
5:  GBP1     Res_immune  9.04 Cytokine
6:  GBP1 Non_res_immune  9.43 Cytokine
7:  CDK6     Res_immune 17.60     IL17
8:   IL6 Non_res_immune 14.50     IL17
9: CXCL8     Res_immune  1.25     IL17

Then, use facet_grid function of ggplot2:

library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_grid(Group~., scales = "free")

Alternatively, you can use facet_grid (or facet_wrap) to have your plot in a single row:

library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_wrap(.~Group, scales = "free")

Reproducible data example

structure(list(gene = c("CCT5", "CCT5", "CDK1", "CXCL5", "GBP1",
"GBP1"), variable = c("Res_immune", "Non_res_immune", "Res_immune",
"Non_res_immune", "Res_immune", "Non_res_immune"), value = c(10.6,
10.5, 9.25, 9.82, 9.04, 9.43), Group = c("Cytokine", "Cytokine",
"Cytokine", "Cytokine", "Cytokine", "Cytokine")), row.names = c(NA,
-6L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

and

structure(list(gene = c("CDK6", "IL6", "CXCL8"), variable = c("Res_immune",
"Non_res_immune", "Res_immune"), value = c(17.6, 14.5, 1.25),
    Group = c("IL17", "IL17", "IL17")), row.names = c(NA, -3L
), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)


Find covariation of different DNA-binding protein binding between two conditions

Answer from @stupidwolf converted from comment:
I would try something like this, for every gene you have a logFC of RNAPII in control versus treatment. Do this for other 2 proteins, and you correlate the logFC. In very rough terms, correlation is the scaled form of covariance, so this will work.
If you have 1000 genes, do correlation between 1000 logFC from RNAPII against 1000 logFC from other proteins.
Answer from @m converted from comment:
You use a covariation stat and then perform Pearson's correlation to standardise the response between -1 to 1. An answer of >0.7 is considered robust and >0.8 highly robust. ITs an interesting question and you are right covariation is the way to answerr it


How to find all WGS assemblies accessions of a species

There seem to be WGS assemblies that can't be found in NCBI Assembly database, e.g.: https://www.ncbi.nlm.nih.gov/nuccore/1779902990. I guess that such assemblies also cannot be found in The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
My current best guess is that each WGS assembly has a "WGS master record" in NCBI Nuccore database. To find all WGS assemblies of a taxon whose uid in NCBI Taxonomy database is 1337, search NCBI Nuccore database using Entrez and the search term "wgs master"[Properties] AND txid1337[orgn:exp].


Error with BWA Mem input having multiple fastq files using cat and process substitution

You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment.

From bwa docs: If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair. BWA will not search through FASTQ files trying to find matching mate IDs.

You will notice that your initial error message is reporting two different read IDs that are improperly paired. Presumably the mispaired read corresponds to the first read in the out-of-order file.


How do I tell featureCounts to ignore fragments that overlap with more than one feature?

If a fragment overlaps multiple genomic features then it's unassignable, which is the appropriate way to handle ambiguous mappings when the resulting counts are then given to tools like DESeq2.


This is really tough to do with highly heterozygous animals. What are your genome assembly stats? Specifically, what is your number of contigs, scaffolds, assembly size, and the N50?

If you have proximity ligation data it will be easier to determine if potential paralogs are truly from different regions of the genome or are just from homologous locations on sister chromatids that ended up in uncollapsed in your final assembly.

If I were you and had a good genome, I would start by annotating gene models using RNAseq data, then blasting each gene in the model against all genes in the model. This would give you potential paralogs to start looking at more closely.


How do I write tests for a snakemake pipeline?

You can run CI on a snakemake workflow when provoding a minimal example dataset. Maybe an example from the snakemake-workflows repo will be helpful: https://github.com/snakemake-workflows/dna-seq-varlociraptor


Using ONT MinION, why is that we cannot get a full length DNA read?

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


The Refseq team and also the NCBI resource coordinators team publish a new paper every few years, so check out the many papers (e.g. here or here), but to answer your 2nd question, non-redundancy here is (I think) defined very strictly as proteins that are identical in terms of sequence and length, so the clustering is trivial, without the need for a sophisticated clustering algorithm as required to detect more remote homologs.


Pipeline for extracting gene from multiple genomes for use in HyPhy selection analyses?

Okay well, because you're using HYPHY, you're going to need a gene tree of every gene you intend to analyze. This requires you to understand the homology relationships between the genomes at hand. Depending on your project, you might already have a list of homologs you're trying to analyze, or you can start from scratch.

If starting from scratch, you will need the proteomes of the species you're interested in (that is, the translation of the coding sequences of each gene, or to put it another way, the concatenated exons). Genome releases should typically have a proteome available for download, along with its corresponding coding sequence.

Next, you need some way to get the homolog groups for the genes in your genomes. You can use a pipeline like OrthoFinder, or iterative tree building and trimming, to infer homolog clusters (OrthoFinder calls them orthogroups but they're the same thing - homologous genes for the species of interest including both single copy orthologs and paralogs).

For this step, if you already have a set of genes and you know the copy numbers, or know that they're single copy orthologs, then it is more easy - you can just pull the sequence from each proteome and go from there.

Once you have your homologs, you need to align the amino acids using a tool like MAFFT, and then align the coding sequences to the amino acid alignment to attain a coding sequence alignment. For this you can use a tool like pal2nal, or check out pxaa2cdn from the phyx package written by some of my collaborators.

If the sequences aren't too divergent, and intron-exon structure is conserved, you could possibly also include alignments of the introns to provide more phylogenetic information. If you do use this approach however, you'd need to align the exons and introns separately. Also, the introns would not be included in the selection analysis.

Next, I would typically reinfer the gene tree of each gene using nucleotide models, which depending on the time, level of selective constraint, etc. can give more signal for the phylogenetic relationships. For this you can use a tree inference program like RAxML. With this tree and your coding sequence alignment, you'll be ready to run analyses in HYPHY (aBS-REL, MEME, FUBAR, etc., depending on your requirements).

Note that if you want, you can also do this in a Bayesian framework to integrate over uncertainty in the model parameters, tree topology, etc. Several Bayesian phylogenetic inference tools (e.g. MrBayes, PhyloBayes), allow you to infer putatively positively selected sites at the same time as tree inference.


Look up a tutorial on DESeq2.  To keep things simple for your first time, I'd reduce your phenotypes to "high" and "low" bins instead of trying to use the cell count numbers as they are.  Then once you understand what the software is doing, you can get fancier.


How should the SAM MD tag match the CIGAR string?

The MD string doesn't apply to soft or hard clipped regions,. so your example read becomes 43G42. Since variant calling using something simplistic like this is only ever going to use the aligned portion, I guess there was never much reason to bother with adding in soft-clipping information. This also means that you need to parse the CIGAR string if you're using the MD string.


How do I increase the minimum dot size in Seurat's DotPlot function?

Did you try to use DotPlot(..., scale.by = "size")?

In contrast to the default scale.by= "radius", this will link the area (==2*pi*r^2), not the radius, of the circles to the fraction of cells expressing the feature.
This corresponds much better to our perception of size and will make differences in low values easier to see.

Reading ?Seurat::DotPlot the scale.min parameter looked promising but looking at the code it seems to censor the data as well.

Since Seurat's plotting functionality is based on ggplot2 you can also adjust the color scale by simply adding scale_fill_viridis() etc. to the returned plot. This might also work for size. Try something like:

DotPlot(...) + scale_size(range = c(5, 10)) # will like warn about supplying the same scale twice

Apart from this, Seurat's plotting system is not very hackable and I find it much easier to extract the relevant data and plot them myself with ggplot2.

Yet another comment: Your plot with the strong differences looks much more convincing to me wrt. to the marker property of these genese than thee cited plot.


Base abbreviations other than ACTGU in sequence file

I usually convert SRA files to FASTQ before doing any further processing, but the tooling for converting SRA files is cumbersome.
To simplify things, I recommend you download the FASTQ file directly. The ENA currently has a nicer interface for doing these things. Here is a link to download the FASTQ files for the SRA sequencing run associated with the SRA experiment you are describing: https://www.ebi.ac.uk/ena/browser/view/SRR11588968
A description of the FASTQ format can be found here: https://en.wikipedia.org/wiki/FASTQ_format


Secretome and membrane receptor profiles

I don't think there is a database with expression data for your particular cell type. I am afraid your "headache" approach doesn't sound so bad...

So get expression data e.g., from GEO for your cell type, and define which genes/proteins you consider as "expressed" (e.g., minimal RPKM value).

Then, you can use GO annotation to find your genes involved with membrane or secretion (or what ever). You can use biomaRt in R for this.

For example, you want to know all "membrane" genes.

> library(biomaRt)
>
> ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
>
> membrane.genes <- getBM(attributes=c('hgnc_symbol',
> 'ensembl_transcript_id', 'go_id'),
> filters = 'go', values = 'GO:0016020', mart = ensembl)

And now you want to know it only from genes that are expressed in your ganglion cells (with dummy example of two genes).

> gene.list <- c("KIR2DL3", "GPR1")
>
> membrane.genes.list <- membrane.genes[membrane.genes$hgnc_symbol ==
> gene.list,]
>
> membrane.only <- membrane.genes.list[membrane.genes.list$go_id ==
> 'GO:0016020',]


filter out FASTQ reads which are shorter

You need to use variable $seq. For example:

bioawk -cfastx 'length($seq)>=259{print "@"$name"\n"$seq"\n+\n"$qual}' test.fq.gz

There are also more convenient tools like seqtk and seqkit. With those, you may

seqtk seq -L260 test.fq.gz


A reverse image search on google leads to the exact code used to produce the second image:

[Xcoarse, Ycoarse] = meshgrid([0 1 2 3], [0 1 2 3]);
[Xfine, Yfine] = meshgrid(linspace(0,3,3000), linspace(0,3,3000));
DataCoarse = [ 1 2 4 1; ...
           6 3 5 2; ...
           4 2 1 5; ...
           5 4 2 3];
DataBicubicFine  = interp2(Xcoarse, Ycoarse, DataCoarse, Xfine, Yfine, 'bicubic');
figure
surf(Xfine, Yfine, DataBicubicFine); shading flat; colormap(jet); view(0, 90)
 hold on; plot3(Xcoarse, Ycoarse, 10*ones(size(Xcoarse)), 'k.', 'MarkerSize', 20)
colorbar
print -dpng -r300 BicubicInterpolationExample.png


Why does cutadapt remove low quality bases from the ends of reads only?

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).


ncbi...eutils/esearch.fcgi: How to search now in the snp database?

eUtilities for dbSNP is in the process of being retired. See this post from NCBI. You may want to check out the new API.
From the NCBI Insights blog post:

Entrez search is available only for human. We will retire all Eutils report formats (RS docsum XML and ASN.1, FASTA, and FLAT text) later this year. Only eSummary XML report will be supported.
Batch Query Service for all organisms will be retired October 2018. Users of human dbSNP data should transition to the new API.
We will place non-human FTP data in the archive directory. Please see the phase-out plan announced last year.


How to create a graph of multi column gene dataset?

As far as I understand, it seems that the first table you provided contains biological/experimental information related to your genes, i.e. the nodes of the network you will build. The second table is the one you need to build the network: each line is an edge. You can use it as input and build the network using the IGraph library which works for both Python (as you requested) or R (as you suggested as an alternative).

Most probably, but I can't say for sure, the first table was just a dataset from which you (or someone for you) recovered a list of nodes that where then submitted to a database, i.e. String, commonly used to recover information related to the interactions between your nodes, i.e. genes. Once you have this info (in the second table) you can build the network!

edit:

if you're not familiar with IGraph, you can have a look at this question describing how to obtain an IGraph graph from a csv file or this other one which uses a different csv file format to build a graph


You can use contrast = c(0.5, 0.5, -0.5, -0.5), which means "the average of the first two coefficients minus the average of the second two.


Python equivalent to R GRanges and IRanges

I am not sure if your question has anything to do with IRanges. If I get it correct, the limitation you have with pyranges now is adding extra columns to the meta data, like adding values to @values slot in granges.
You can do it using insert or setattr :
import pyranges
import pandas as pd
import numpy as np
gr = pyranges.random(10)

setattr(gr,'alt',np.random.choice(['A','T','G','C'],10))

gr.insert(pd.DataFrame({'ref':np.random.choice(['A','T','G','C'],10)}))

+--------------+-----------+-----------+--------------+------------+------------+
| Chromosome   | Start     | End       | Strand       | alt        | ref        |
| (category)   | (int32)   | (int32)   | (category)   | (object)   | (object)   |
|--------------+-----------+-----------+--------------+------------+------------|
| chr1         | 27431195  | 27431295  | -            | C          | G          |
| chr2         | 197045893 | 197045993 | -            | A          | G          |
| chr4         | 86316012  | 86316112  | -            | C          | A          |
| chr8         | 4560598   | 4560698   | -            | G          | G          |
| ...          | ...       | ...       | ...          | ...        | ...        |
| chr11        | 130395526 | 130395626 | -            | A          | A          |
| chr12        | 28557156  | 28557256  | +            | T          | G          |
| chr13        | 55519337  | 55519437  | -            | T          | C          |
| chr14        | 1807028   | 1807128   | +            | G          | A          |
+--------------+-----------+-----------+--------------+------------+------------+
Stranded PyRanges object has 10 rows and 6 columns from 10 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.

Or directly adding them like this (thanks to @ChrisRands for pointing it out):
gr.alt = np.random.choice(['A','T','G','C'],10)


Classifying samples based on marker gene expression

I would consider using gene expression signatures to classify samples (especially cancer subtypes but the same principles apply to other problems of this type) one of the classic problems of bioinformatics. Quite a lot of work has been done on methods to derive gene sets that provide good classification performance. This is slightly different from your problem since you already have a gene signature but it may still prove useful.

These methods will typically fit a model that selects a (small) number of genes from genome-wide expression data that distinguish between the cell types/conditions in question, i.e. they derive a gene signature. The resulting model then allows for the classification of new samples. I've had success using GeneRave for this purpose (but note that this was designed for microarray data, I haven't used it with RNA-seq data and don't know how well it holds up there). A more recent paper relating to this issue can be found here.

So how does that help you? One option would be to fit one of these classifiers to gene expression data for the genes you already know to obtain a model that can then be applied to new samples automatically.


Update

As pointed out by neilfws:

According to Wikipedia and Practical bioinformatics by Michael Agostino: The name "Entrez" means "Come in!" in French (pronounced “on-tray”); it was chosen to reflect the spirit of welcoming the public to search the content.

Old (probably incorrect) theory

The first version of Entrez database was distributed by NCBI in 1991 (on CD-ROM). At that time, it consisted of nucleotide sequences, protein sequences and associated citations and abstracts from MEDLINE.

So even at the early stage it contained objects (entries) of very different types (sequences and abstracts). It is hard to find a word which will accurately and briefly describe all the types of data.

While I don't know the correct answer and this is purely a speculation, my guess is that Entrez is an intentionally misspelled word "entries"; which accurately describes all the datatypes in the database.

Using a misspelled word as a name is not uncommon, e.g. Google comes from misspelled "googol" and Flickr is a misspelled version of "flicker".


Interpreting quantitive outputs from maximum likelihood phylogenetic trees

In summary,

RAxML_bipartitions.output_bootstrap.tre

Is the only file of interest. The reason this is true in this context is really complicated and you have to understand the statistics of likelihood and how they are interpreted within phylogeny to understand why. This file is simply the final output of a non-parametric bootstrap analysis performed by maximum likelihood.

What on earth is a non-parametric boostrap?

A non-parametric bootstrap is resampling each alignment position with replacement. Thus if we have alignment positions 1,2,3,4,5
A bootstrap resample for 2 replicates might be,

Replicate 1

1,1,3,5,2

Replicate 2

4,2,5,2,1

The ML algorithm will make trees of replicates 1 and 2 and find the consensus between them. If you think about it in any other context a bootstrap replicate is pretty meaningless because it no longer reflects the true biological sequence. Thus information on how the consensus was derived, are not really of interest to us providing we are confident this has been done correctly, viz. RAxML_bipartitionsBranchLabels.output_bootstrap.tre and RAxML_bipartitionsBranchLabels.output_bootstrap.tre

So why is this output of limited use?

There are situations to some investigators this information is useful, but assess the robustness of a tree topology its not needed. The only thing we want is a phylogram (bestTree) with the bootstrap values superimposed on them. We really don't need complicated stuff such as the tree to be represented for example as a polytomy (non-bifurcating tree) because we can just read the bootstraps to make that deduction (values >> 75%). In addition, there is not perfect consensus what boostrap value constitutes robustness, but generally most agree >80% is robust.

What output files have useful information in them?

The information that is important are the files associated with "bestTree", that was the single maximum likelihood tree performed on the intact native sequence. The "info" file for this contains 3 really important parameters:

-lnL ... very important!!
Gamma distribution parameter "alpha",
PINVAR, proportion of invariant sites,

-lnL is the highest log-likelihood (probability) of the phylogeny. It is usually a very small number for which where is an enormous amount of theory over it.

Alpha parameter of the gamma distribution this is the shape parameter of the mutation rate, if it is very low (<1) the distribution of mutations across the alignment is very tight clustered and approximates to a negative binomial distribution. This means some sites don't mutate at all and a small number of sites mutate alot. If it is very large >200 (which is never observed) it approximates to the Poisson distribution, meaning the mutation distribution is randomised across the alignment.

PINVAR this is a straight percentage/frequency and simply means the sites that don't mutate.

How are they calculated?

PINVAR and alpha are not emperically calculated, i.e. if you look at an alignment and say 'no mutations at that position', PINVAR would of course agree but may consider other invariant depending on the phylogeny. These parameters are calculated by maximum likelihood and you can begin to see why the calculation takes so long ... alpha and PINVAR affect the tree topology (which affects -lnL), but the topology affects alpha and PINVAR. Thus, is a multidimensional search of tree and parameter space.

So what stuff do I report in my Results?

Anyway reporting -lnL is good technique and shows the reader you've done maximum likelihood, citing PINVAR and alpha from gamma distribution helps ('Methods' parameters were calculated reiteratively under maximum likelihood). This is only useful for bestTree. The -lnL, PINVAR and gamma's alpha are also calculated for every single bootstrap replicate, but these values are of limited to use, because we have resampled the data, only the consensus tree counts... Obviously presenting the bootstrapped phylogram is extremely important.

Welcome to the technical world of phylogeny!

The amino acid matrix you used BTW .. LG is in vogue right now.

How do I do it?

When I do this stuff its via Biopython and ETE3, I capture the values within the pipeline and don't examine the output files of RAxML because I generate my own.


different results coming from biomart online and biomaRt R library

Sorry, this was my mistake in the last question. To search down the ontology, rather than just for the specific association with a term, biomaRt needs a different filter: 'go_parent_term'.

Try:

gene <- getBM(attributes = c('external_gene_name'), filters = 'go_parent_term', values = 'GO:0030098', mart = ensembl)


In calculating the retention index, why do we use the character state with the lowest frequency?

The maximum number of steps (or changes) is the number of taxa with state 1 or 0, whichever is smaller. Conceptually, this is the number of character changes if each taxon evolved its state independently of the other taxa.

If for a given character the character state 0 occurred twice and character state 1 occurred five times, 2 would be the maximum number of changes.


DESeq2 complicated design - effect of replicated samples

I agree with @devon-ryan that the design matrix looks reasonable, i.e. that it's okay to treat the new samples as biological replicates of the original batch, and that the batch number should be included in that matrix (assuming DESeq2 doesn't complain about it) to filter out most of the batch effects.

DESeq2 tends to be conservative in its differential expression calculations. While it's true that more replicates in one sample will increase the number of differentially-expressed genes that are statistically significant (because the biological variation can be better modelled), I wouldn't expect this to change the overall picture.

However, whole-genome differential expression testing should not be used on its own to demonstrate biological relevance. I recommend as a first-pass sanity check that people compare any differential expression results to what is seen on a genome browser to identify any quirky things going on (e.g. expression across only a part of the gene). Ideally, results should be confirmed with targeted analyses on the same (or similar) biological samples, either at the transcript level (e.g. qPCR) or at the protein level (e.g. Western Blot). When aggregating evidence from multiple analyses, the specific way in which one of those analyses was done is less important.

Aside: It'd obviously be better [in the future] to make sure that your sample replicates are both robust to drop-out, and allow for outliers to be noticed, reducing the batch-effect problem. In practise this means a minimum of four replicates per informative condition, but because funders / researchers are well known for only choosing the minimum, I'd just state that six replicates per sample as what is needed to produce good results.


How to proceed after contigs.fa file obtained from VELVET assembly

The contigs.fa file is a FASTA formatted file of the final contigs generated by velvet.  I would start assessing the quality of the assembly that was generated. Here are some questions to get you started.  Really the next steps depend on how you plan to use the assembly.

What is the distribution of contig lengths? Average and N50 contig length

How does the assembly compare to the reference?

How much of the reference does the assembly cover?

When you align the original reads back to the contigs how well are they covered?

Do any of the read alignments not make sense?
Are there sections of the contigs that are not covered?


Extracting all information about a sample when using xtract from e-utilities

You can give pysradb a shot.
study_accession experiment_accession experiment_title                                                   experiment_desc
                                      organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument           total_spots total_size  run_accession run_total_spots run_total_bases
SRP178555       SRX5243234             GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq    GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245909                    Illumina HiSeq 2500  19052916     446793895  SRR8435699    19052916         971698716
SRP178555       SRX5243233             GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq    GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245908                    Illumina HiSeq 2500  18156224     428457264  SRR8435698    18156224         925967424
.. [truncated]


I like the idea but I’m not aware of a strong Lisp presence anywhere in bioinformatics. I think a good starting point would be the paper How the strengths of Lisp-family languages facilitate building complex and flexible bioinformatics applications by Khomtchouk &al (and the Hacker News discussion).

Excluding Clojure from consideration will unfortunately close off a great many modern resources (although I fundamentally agree with your technical reason for doing so). Apart from this, Scheme (in particular R7RS) is probably the most-used Lisp dialect, ahead of Common Lisp. Or at least the most actively developed one. As Pierre mentioned, there’s also GNU Guile but, to be honest, despite it being the default for GNU software tools I don’t think it’s frequently used outside GDB.

That said, Clojure’s startup time is not anywhere near 10 seconds on my system. It’s not exactly instantaneous but not much worse than other dynamic languages, despite being backed by the rather sluggish JVM.


There seems to be a dividing line around CoV=0.9.  Perhaps this is a good place to start for a cutoff.


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

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


DiffBind: Model must include count data for contrasts.

The problem was that I was using wrong files. I ran a pipeline from here and it generates several folders with many files. The .bam files that should be used are located in the aligned_mm10 folder and are with the extension: .pe.q10.sort.rmdup.bam.


Link Prediction in Bipartite Networks between biological pathways and drugs

I am not sure if I understood correctly how you classify  the drugs. But what you attempt to do is similar to what they do here but in that article they use the drugs targeted to a molecule/pathway to find combination of drugs that are better for the disease (not an alternative drug), see the first image:

In that article they mention the pathway cross-talk inhibition (PCI), you might want to calculate it for your pathways to assess the relative reduction of the pathway efficiency.


Scatter tasks around predefined batches inside a WDL workflow

Here's an alternative to using a custom task, albeit still rather unappetizing:

workflow groupby {
  input {
    File? nil
  }

  Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
  Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

  scatter (y in ["A", "B"]) {
    scatter (p in zip(X,Y)) {
      File? x_or_nil = if (p.right == y) then p.left else nil
    }
    Pair[String,Array[File]] group = (y, select_all(x_or_nil))
  }

  output {
    Array[Pair[String,Array[File]]] ans = group
  }
}

(It'll be slightly cleaner with the upcoming None literal instead of dummy nil)

FWIW, here's a lengthy discussion about the idea of allowing developers to define custom library functions which just run a small task that manipulates data structures (hinting that it needn't necessarily offload to a worker node, with the huge latency that typically entails).


Its easy,

Download the data from NCBI in fasta format (its an option in one of the pull down headers).
Load it into the latest version of MEGA, version 7(?), perform the alignment therein.
Export from MEGA as Nexus format and upload it into Beast

You could of course perform the alignment in Clustalx, Clustalo or Muscle or T-Coffee, (hmmm I could write an alignment program just to name it T-4two (sorry British joke)). If the output option of "nexus" was not present in the alignment program of your choice, then again upload the alignment into MEGA and export it as nexus.

Anyway, I could go into specific detail about the various points and clicks therein, but I think its easy to figure out. Most phylogeny programs will do an automatic fasta->nexus conversion, but anyway.

Its easy, Beast is not trivial to use properly however and ensuring convergence is not at all easy, in addition the multiple simulation a priori models, whilst powerful, give different answers. To try and explain this a different way, it is easy to get an answer in Beast, but difficult to get a robust answer from it and thats the bit I would give serious consideration about. It does amaze how widespread the program is (I guess its the nice GUI), given how complicated it is to obtain robust result.

Just to add you could try the online tool here for the fasta to nexus conversion:  . I've not tried it, its just an after thought. If you do Python, Biopython will do it, but its just as easy to code.


Simplest way to work out structural variant type?

GRIDSS's author had posted his own solution on Github. The code is:

simpleEventType <- function(gr) {
    return(ifelse(seqnames(gr) != seqnames(partner(gr)), "ITX", # inter-chromosomosal
           ifelse(gr$insLen >= abs(gr$svLen) * 0.7, "INS",
           ifelse(strand(gr) == strand(partner(gr)), "INV",
           ifelse(xor(start(gr) < start(partner(gr)), strand(gr) == "-"), "DEL", "DUP")))))
}


Difference between isoforms and paralogs in transcriptoms?

If a reference genome assembly exists, then reference-guided assembly with Trinity will keep paralogs, because it treats each genomic region as a different assembly set.
If a reference genome doesn't exist, paralogs will probably result in the construction of non-representative "isoforms" containing variants present in multiple genes. Inchworm does a greedy assembly of reads, so if a chimeric transcript can be constructed from the reads (which is likely with paralogous genes), then it might end up being fed to the subsequent steps. Whether that is actually emitted at the end by Butterfly will mostly depend on how the paths through the putative assemblies are traversed, which probably depends on the tolerance of Bowtie to errors.
If paralogous transcripts are sufficiently similar that the differences are within the error tolerance for assembly, then I would expect that they'd be collapsed into a single transcript.


How to convert a list of genes to one element per row?

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

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


The PCHeatmap function (renamed DimHeatmap in Seurat v3) can be used to help determine the number of principal components to use in downstream analysis, as well as to visualize the top genes contributing to each PC.

Both cells and genes are ordered by their PC scores, and by default the 15 genes with the highest and 15 genes with the lowest PC loadings are displayed for each PC.

In your case with PC1 and PC2, my conclusion would be that PC1 and PC2 both strongly separate different populations of cells. For PC1 there is a smaller population with very strong separation, and PC2 splits the cells roughly in half.

You can also set the cells.use parameter (renamed cells in Seurat v3) to only show the top n cells for each PC which can help with visualization.


europmc and tidypmc R library for extracting or making metadata from publication

The idea for this code is to first convert PIDs to PMCIDs, then run the tidypmc in a loop over the PMCIDs. The only problem is that tidypmc failed to retrieve tables from most of the IDs in your example list.
library(tidyverse)
library(tidypmc)
library(httr)
library(jsonlite)
example_pids <- c(30024784, 30555165, 30510081, 31688884, 31516032, 28588019, 29286103) %>% as.character()

#-- Convert to PMC ids
convertPIDtoPMCID <- function(pids) {
    #-------- Make API request
    pids4query <- paste(pids, collapse = "%0D%0A")
    idconv_req <- paste0("https://www.ncbi.nlm.nih.gov/pmc/utils/idconv/v1.0/?ids=", pids4query, "&idtype=pmid&format=json&versions=no&showaiid=no&tool=&email=&.submit=Submit")
    pids_json <- GET(idconv_req) %>% content("text") %>% fromJSON()
    #-------- Get info from JSON
    pid2pmc <- pids_json$records %>%
        select(pmcid, pmid) %>%
        as.data.frame()
    rownames(pid2pmc) <- pid2pmc$pmid
    pmcids <- pmcids <- pid2pmc[pids, "pmcid"]
    return(pmcids)
}
example_pmcids <- convertPIDtoPMCID(example_pids)

#-- Try to get data with tidypmc
pub_tables <- lapply(example_pmcids, function(pmc_id) {
    message("-- Trying ", pmc_id, "...")
    doc <- tryCatch(pmc_xml(pmc_id),
                    error = function(e) {
                        message("------ Failed to recover PMCID")
                        return(NULL)
                        })
    if(!is.null(doc)) {
        #-- If succeed, try to get table
        tables <- pmc_table(doc)
        if(!is.null(tables)) {
            #-- If succeed, try to get table name
            table_caps <- pmc_caption(doc) %>%
                filter(tag == "table")
            names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ")
        }
        return(tables)
    } else {
        #-- If fail, return NA
        return(NA)
    }
})
names(pub_tables) <- example_pids

#-- Inspect results
pub_tables$`30555165`$`Table 1 - Patient demographic and baseline characteristics`
pub_tables$`29286103`$`Table I - Sample summary.`

Tables will require quite a bit of tidying after this. Good luck!


Getting species from ensembl sequence id

The lookup/id endpoint will get it for you. Of you can just look up the three letter species code, in this case MOD, on the list in the documentation.


How to safely and efficiently convert subset of bam to fastq?

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter list), store the properly oriented sequence/quality/read name in a buffer, and then dumps that buffer entry to disk once the mate is found. This should be reasonably performant (hey, it's python, to don't expect too much), though if you happen to have indexed BAM files then one could think of ways to make this run faster.

Do check the output, since I've only run one test.


Clustered column chart for single cell data

You can make use of table() and prop.table() to calculate the number and proportions of your cell types with respect to tissue of origin.

prop.table() requires a table as input so you can start with something like prop.table(table(object@meta.data$cell_type, object@meta.data$tissue). This will create an objet of class table.

Next step is to create a data frame with as.data.frame(your_table_object) and then to convert it to "long format", for example with tidyr::pivot_longer(). You can feed these data in long format to ggplot().


How to interpret the GROMACS MD Simulation RMSD Graph?

The opposite-ish.

Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.

RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.

Structure RSMD

If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.

However, the command you used was based on gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg from the page you linked. Key here is the index -n index.ndx, which limits it to the ligand. Before this happened gmx trjconv aligned the protein in the trajectory, so that common mistake was not done.

Ligand RSMD

In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (gmx energy) will tell you how strongly does it bind.

The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD.
To explore this further you can:

look at the trajectory
in PyMOL using the command split_states twice, once for the first and then the last state, followed by align of the two
If it were the whole protein, looking at the per bond matrix -bm option in Gromacs rms and spotting the change.

3Å??

I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.

Did the energy minimisation converge properly?
If it is because the proteins themselves wobble substantially, revisit gmx trjconv with a focus on the active site minus the ligand.
Is the ligand binding but leaving one substituent group to wobble?

Energetics

As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k_off and k_on and their ratio K_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k_on.


Are there certain alignment methods/tools which perform better with a high density of indels?

Minimap2 attempts to solve a lot of problems that are not properly dealt with by other mapping software. It's written by the same author as BWA, but designed to deal with situations where the reference and query sequences are quite different (e.g. for noisy long reads).


I have really skewed RNA-seq data, what's the best way to normalise it? Preferably in R!

Welcome to RNA-seq. This is not at all unusual, as RNA-seq raw data are not normally distributed. People typically approximate it via the Negative Binomial distribution. I suggest you check existing packages as Devon suggests in his comment for normalization and differential expression (or anything else). Typical choices are DESeq2, edgeR or limma, but there are others as well.
This guide might help you get started. Please also use google to explore existing threads both here at SE and over at biostars.org.
Most basic questions have already been discussed there and are informative to get started.


tool for the prediction of impact of mutations on protein localization

Do you have the protein sequence after the mutation available with you?
If so then maybe you can BLAST them and find out where the best hits are localized using GO annotation (cellular components). I am aware that these won't be your exact protein sequence but they can at least give you a rough idea.


GenomicRanges::countOverlaps seems to be what you’re after:

position_range = GRanges(position$chrom, IRanges(position, position, width = 1))
ranges_at_position = countOverlaps(position_ranges, granges)


Detecting portions of human proteins with high degree of microbial similarity

Sounds like precisely the job BLAST was developed for. Now, which flavor will depend on what you want to do and what data you have available. Some options:

PSI-BLAST: this is usually the best choice if you are trying to find protein homologs. It works by building a hidden markov model describing your query sequence and using that model to query a database of proteins. The advantage is that it is run in multiple iterations, giving you the chance to add or remove results (so you add the ones that are true positives and remove false ones), eventually building a pretty good model of your protein. This is far moer powerful than a simple homology-based approach since proteins work via protein domains and simple homology is not as important as specific conserved functional residues.

For this, go to the NCBI protein blast page and select PSI-BLAST:

BLASTp: Simple protein-protein blast. It will identify homologous proteins based on sequence similarity. Whether or not that also implies functional homology is not that simple and will depend on each case you investigate.

As above, go to the NCBI protein blast page, but this time use the defaults.
tBLASTn: this is a tool that takes protein sequences as input and compares them to a database of DNA which is dynamically translated in all 6 possible reading frames. Very good for finding homologous sequences when you don't have well annotated protein information for the target species. It has the benefit of being more sensitive and able to find more distant homologies than basic nucleotide BLASTn and the go-to approach when your target species is distant and not well annotated.

NCBI's tBLASTn page.

All of these can be run online through the NCBI's BLAST page. If you want to investigate hundreds of proteins, I suggest you install blast locally. You can then either download the relevant target sequences from NCBI and rebuild the blast database locally (if so, I suggest you ask a new question about how to do that) or, use NCBI's remote blast client which lets you use a locally stored query file and will run blast on the NCBI's servers.

Now, these programs will return what are known as High Scoring Pairs (HSPs), the regions of your query sequence(s) that align well to the target. There are various options you can play with to improve the sensitivity or the specificity but a discussion of those would require far more details about what you're doing and would be best in a new question as well.

Once you have your HSPs, you can relatively easily parse them to select regions with a given range of sequence similarity values and of a specific length. Once again, that would be better discussed in a separate question once you have your results and can show an example.


Quast duplication ratio and mismatches percent

I took a quick look at the quast code that I believe is responsible here:
def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
    indels_info = IndelsInfo()
    genome_mapping = {}
    for chr_name, chr_len in reference_chromosomes.items():
        genome_mapping[chr_name] = [0] * (chr_len + 1)
    with open(used_snps_fpath, 'w') as used_snps_f:
        for chr_name, aligns in ref_aligns.items():
            for align in aligns:
                ref_pos, ctg_pos = align.s1, align.s2
                strand_direction = 1 if align.s2 < align.e2 else -1
                for op in parse_cs_tag(align.cigar):
                    if op.startswith(':'):
                        n_bases = int(op[1:])
                    else:
                        n_bases = len(op) - 1
                    if op.startswith('*'):
                        ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                        if ctg_nucl != 'N' and ref_nucl != 'N':
                            # MISMATCHES INCREMENTED HERE
                            indels_info.mismatches += 1
                            if qconfig.show_snps:
                                used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += 1
                        ctg_pos += 1 * strand_direction
                    elif op.startswith('+'):
                        # INDELS INCREMENTED HERE
                        indels_info.indels_list.append(n_bases)
                        indels_info.insertions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = '.', op[1:].upper()
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ctg_pos += n_bases * strand_direction
                    elif op.startswith('-'):
                        indels_info.indels_list.append(n_bases)
                        indels_info.deletions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = op[1:].upper(), '.'
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += n_bases
                    else:
                        ref_pos += n_bases
                        ctg_pos += n_bases * strand_direction
                if align.s1 < align.e1:
                    for pos in range(align.s1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
                else:
                    for pos in range(align.s1, len(genome_mapping[align.ref])):
                        genome_mapping[align.ref][pos] = 1
                    for pos in range(1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
            for i in ns_by_chromosomes[align.ref]:
                genome_mapping[align.ref][i] = 0

    covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
    return covered_ref_bases, indels_info

This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.
This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).
From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast used_snps file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.


What is the difference between FASTA, FASTQ, and SAM file formats?

Let’s start with what they have in common: All three formats store

sequence data, and
sequence metadata.

Furthermore, all three formats are text-based.

However, beyond that all three formats are different and serve different purposes.

Let’s start with the simplest format:

FASTA

FASTA stores a variable number of sequence records, and for each record it stores the sequence itself, and a sequence ID. Each record starts with a header line whose first character is >, followed by the sequence ID. The next lines of a record contain the actual sequence.

The Wikipedia artice gives several examples for peptide sequences, but since FASTQ and SAM are used exclusively (?) for nucleotide sequences, here’s a nucleotide example:

>Mus_musculus_tRNA-Ala-AGC-1-1 (chr13.trna34-AlaAGC)
GGGGGTGTAGCTCAGTGGTAGAGCGCGTGCTTAGCATGCACGAGGcCCTGGGTTCGATCC
CCAGCACCTCCA
>Mus_musculus_tRNA-Ala-AGC-10-1 (chr13.trna457-AlaAGC)
GGGGGATTAGCTCAAATGGTAGAGCGCTCGCTTAGCATGCAAGAGGtAGTGGGATCGATG
CCCACATCCTCCA

The ID can be in any arbitrary format, although several conventions exist.

In the context of nucleotide sequences, FASTA is mostly used to store reference data; that is, data extracted from a curated database; the above is adapted from GtRNAdb (a database of tRNA sequences).

FASTQ

FASTQ was conceived to solve a specific problem arising during sequencing: Due to how different sequencing technologies work, the confidence in each base call (that is, the estimated probability of having correctly identified a given nucleotide) varies. This is expressed in the Phred quality score. FASTA had no standardised way of encoding this. By contrast, a FASTQ record contains a sequence of quality scores for each nucleotide.

A FASTQ record has the following format:

A line starting with @, containing the sequence ID.
One or more lines that contain the sequence.
A new line starting with the character +, and being either empty or repeating the sequence ID.
One or more lines that contain the quality scores.

Here’s an example of a FASTQ file with two records:

@071112_SLXA-EAS1_s_7:5:1:817:345
GGGTGATGGCCGCTGCCGATGGCGTC
AAATCCCACC
+
IIIIIIIIIIIIIIIIIIIIIIIIII
IIII9IG9IC
@071112_SLXA-EAS1_s_7:5:1:801:338
GTTCAGGGATACGACGTTTGTATTTTAAGAATCTGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII6IBI

FASTQ files are mostly used to store short-read data from high-throughput sequencing experiments. The sequence and quality scores are usually put into a single line each, and indeed many tools assume that each record in a FASTQ file is exactly four lines long, even though this isn’t guaranteed.

As with FASTA, the format of the sequence ID isn’t standardised, but different producers of FASTQ use fixed notations that follow strict conventions.

SAM

SAM files are so complex that a complete description [PDF] takes 15 pages. So here’s the short version.

The original purpose of SAM files is to store mapping information for sequences from high-throughput sequencing. As a consequence, a SAM record needs to store more than just the sequence and its quality, it also needs to store information about where and how a sequence maps into the reference.

Unlike the previous formats, SAM is tab-based, and each record, consisting of either 11 or 12 fields, fills exactly one line. Here’s an example (tabs replaced by fixed-width spacing):

r001  99  chr1  7 30  17M         =  37  39  TTAGATAAAGGATACTG   IIIIIIIIIIIIIIIII
r002  0   chrX  9 30  3S6M1P1I4M  *  0   0   AAAAGATAAGGATA      IIIIIIIIII6IBI    NM:i:1

For a description of the individual fields, refer to the documentation. The relevant bit is this: SAM can express exactly the same information as FASTQ, plus, as mentioned, the mapping information. However, SAM is also used to store read data without mapping information.

In addition to sequence records, SAM files can also contain a header, which stores information about the reference that the sequences were mapped to, and the tool used to create the SAM file. Header information precede the sequence records, and consist of lines starting with @.

SAM itself is almost never used as a storage format; instead, files are stored in BAM format, which is a compact binary representation of SAM. It stores the same information, just more efficiently and, in conjunction with a search index, allows fast retrieval of individual records from the middle of the file (= fast random access). BAM files are also much more compact than compressed FASTQ or FASTA files.

The above implies a hierarchy in what the formats can store: FASTA ⊂ FASTQ ⊂ SAM.

In a typical high-throughput analysis workflow, you will encounter all three file types:

FASTA to store the reference genome/transcriptome that the sequence fragments will be mapped to.
FASTQ to store the sequence fragments before mapping.
SAM/BAM to store the sequence fragments after mapping.


Find a cutoff value for genes that are expressed in single cell RNA-seq?

Not sure if this has anything to do with your question but I think the way how Monocle approach this question sort of makes sense to me. The normalization methods for Seurat and Monocle are somehow similar (correct if I am wrong, this is simply my personal experience), as @kohlkopf mentioned.

For monocle specifically, it has this function where you can display dispersion v.s. expression of the genes to decide the cut off. Here is an example of this plot that I generated when processing the PBMC data downloaded from the Seurat tutorial.

You can assign the cut-off manually but I doubt I answered your question on how to decide the cut-off value. Maybe you need to try a set the expression cut-off from a range to experiment what value may potentially answer a biologically relevant question in your context.

If you are interested in how monocle approach this, I think the Monocle introduction and Dave Tang's Blog are excellent in explain how it works.

Again, you are right that 0 could be sequencing error. But I would also suggest that, if you are looking at genes at very small expression level, scRNA-seq, at least droplet-based technique, could be really difficult in such context. I happened to hear someone trying to find lncRNA expression talking about how a "0" sometimes may be a true "0" due the low expression nature of lncRNA in scRNA-seq, etc.

Anyhow, not sure if I answer the question but I do hope it helps you somehow. I am not able to make comments but this is the best I can do.


How I deal with this kind of gene expression comparision

No. If the only difference between Tumour_batch1 and Tumour_batch2 was the library prep/sequencing run, then all this will tell you is if the particular sample of RNA taken from the tumour is different from the particular sample of RNA taken from normal. It almost certainly is.

Indeed if you did technical replicates of two RNA samples, both taken from normal, you would almost certainly find differences.

What you want to find is the differences in RNA level between tumour in general and normal in general. To do that you need to make an estimate of the mean and variance of RNA levels in the tumour and in the normal. To do that you need multiple samples from each.

If you need differences within a single patient, it would be more honest just to compute the fold changes without any implication of statistical certainty. You would then need to think about alternate schemes for measuring the confidence in inferences drawn from the data when you consider the actual biological questions.


Produce a single sequential FASTA sequence out of BAM

If you just want the human chromosomes in a FASTA format, why not downloading directly the individual chromosome files (chr*.fa.gz)?

If you don't rely on the "official" releases (for whatever reason), then what you need is first to assemble your reads into genomes, possibly (but not surely) going down to a chromosome-level resolution. Different approaches for assembly are possible/available, mainly depending on (1) what do you biologically need to answer and (2) the type of sequencing data you have - we'll need more details to be able to help you in this direction.

The samtools fasta utility will "just" convert your reads (likely given from a sequencing facility) from an alignment (BAM) format into a reads (FASTA) format, there is no concept of chromosomes here. You won't get assembled chromosomes this way, just extract the sequence in FASTA format from your alignment file.


Can I automate CLUSTALO and output alignment sequence identity?

So, as the first reply said, there's sort of multiple questions here.

For filtering things out from a table, I would absolutely use pandas. https://swcarpentry.github.io/python-novice-gapminder/08-data-frames/ Pandas is a way of managing dataframes in python. It looks like what you want from your table is to remove certain rows, if the content of those rows at a certain column matches certain IDs from a list (pathogen IDs). That's the sort of thing the pandas library is good at.

For managing workflows, I use either bash or snakemake. Snakemake is a python thing. It has some advantages, such as being able to start from an intermediary file (for example, if a workflow is "sequences -> multiple alignment -> phylogenetic tree", it can easily start from the multiple alignment). Also deals with branched workflows quite well. Can be easily installed with conda (anaconda is absolutely a lifesaver for bioinformatics work imo). I still use bash for simple workflows, because it's slightly simpler/faster.

For potential filtering - If you want to keep the top hits as they are, but filter the lower-scoring hits based on whether they have a PDB ID or not, I'd personally look into Entrez queries. They're a way of making your search more specific. Here they're shown for WebBlast, but they work quite well for command line blast too. https://www.youtube.com/watch?v=bxx5uaKjMa8


how can i calculate p value of miRNA raw reads counts using R?

I would start by re-formulating the question in this post: Your end-goal seems to the differential expression analysis (DEA) of your samples (and you have to have to levels for this like treatment-control, mutant-WT, etc ...). P-values are context-dependent and should be calculated based on your research question/hypothesis.

DEA has been addressed by many researchers and R has a lot of packages to help you with this; DESeq2, edgeR, limma, etc ...

I would go to the Biconductor website, search for one of the three above, and go through their rich vignettes and would come back here for specific questions.


Assign cell types to groups of cells based on their gene expression profiles

I think the analysis will involve two steps:

1) Compile a reference list of cell types, with a list of genes expressed in each of them, based on the literature. C. elegans has only ~1000 cells, so I'm sure it's been well researched. Mark which genes are just expressed in that particular cell type, and which are uniquely expressed, this latter will be more useful markers.

Resources to do that:

http://openworm.org/science.html#other
http://www.wormatlas.org/
http://www.wormbase.org
Pubmed.
Your lab's knowledge.

This will involve mostly manual work, but maybe you can find some list.

For example, platelets express VWF. And activated platelets express SELP.

2) Use one of the single-cell RNA-seq packages (e.g. Seurat) to cluster cells and to find the "marker" (DE) genes for each cluster (vs. all other clusters).

For each cluster, look up the marker genes in the reference list and infer the cell type.

You may need to infer the cell type from the functions and descriptions of the genes.
Also see Retrieve detailed gene descriptions

Keep in mind that you will have many confusing doublets in the dataset expressing characteristics (genes) of two cells.

You may also need to re-do t-SNE and clustering on a subset of cells, e.g. blood/immune cells, in an iterative, hierarchical way.

Do you have information on which cells are from the same worm? Once you have a reference list, you can try and identify cells of just one worm.


What is "unmapped read segments" in the output of samtools idxstats?

An unmapped read segment is any entry with bit 4 of the flag field set (i.e., it's unmapped). Assuming the reference to which you aligned reasonably matches the sample from which the reads were generated, less than ~10% of to total reads should be unmapped.
Please note that BAM files have no idea what a "feature" is. They only contain information about reads and the chromosomes/contigs to which they align.


removing nextera transposase adapters, cutadapt

The adapter sequence you have googled is the sequence on the adapter primer. This works when you want to remove primer-dimers.

With transposase adapters or ATAC seq, you have very short fragments and you sequence into the adapter. So what you is sequenced is actually the reverse complement. See below image from trimmomatic manual

So if you look at what worked, and the sequence you googled, they are reverse complement of one another.

You can also see this in the NexteraPE-PE.fa, Trans1 and Trans1_rc @Haci recommended:

more NexteraPE-PE.fa
>PrefixNX/1
AGATGTGTATAAGAGACAG
>PrefixNX/2
AGATGTGTATAAGAGACAG
>Trans1
TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
>Trans1_rc
CTGTCTCTTATACACATCTGACGCTGCCGACGA
>Trans2
GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
>Trans2_rc
CTGTCTCTTATACACATCTCCGAGCCCACGAGAC

To answer your question, if most of your fragments are short, and there's not a lot of primer dimers, then yes you are doing the right thing.


Filter bam using SNP list in bed format with minimum mapping quality and base quality

Thanks Devon ;-)
using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

a bed file 'input.bed'

and the following snippet:

private IntervalTreeMap<Boolean> treeMap=null;
private final int BASE_QUALITY= 20;
@Override
public Object apply(final SAMRecord record) {
if(record.getReadUnmappedFlag() || record.getCigar()==null) return null;
if(treeMap==null) {
    treeMap = new  IntervalTreeMap<Boolean>();
    final String bedStr;
    try {
        bedStr = IOUtil.slurp(new java.io.File("input.bed"));
        }
    catch(java.io.IOException err)
        {
        throw new RuntimeIOException(err);
        }
    Arrays.stream(bedStr.split("[\n]")).
        filter(L->!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L->L.split("\t")).
        map(A->new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I->treeMap.put(I,true));
    }

byte quals[]= record.getBaseQualities();
int refpos= record.getStart();
for(CigarElement ce:record.getCigar()) {
    CigarOperator op= ce.getOperator();
    if(op.consumesReferenceBases())
        {
        if(op.consumesReadBases()) {
            for(int i=0;i< ce.getLength();i++)
                {
                Interval t = new Interval(record.getContig(),refpos+i,refpos+i);
                if(!treeMap.containsOverlapping(t)) continue;
                int readpos = record.getReadPositionAtReferencePosition(refpos+i);
                if(readpos==0 ) continue;
                readpos--;
                if(readpos>=quals.length)  continue;
                if(quals[readpos]>BASE_QUALITY) return true;
                }
            }

        refpos+=ce.getLength();
        }

    }
return null;
}

usage:

samtools view -bu -L input.bed input.bam |\
   java -jar dist/samjdk.jar --body -f snippet.code


Generating multiple NCBI URLs with eutils

I don't know anything about php, but when I copy your code into here it isn't happy about the variable name $query.

This code appears to produce valid urls, but I can't make the get part work:

$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene";$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae"); $links = array();

foreach ($gene_symbols as $key => $GS){
   $SP = $species[$key];
   $SearchTerm = "$GS AND $SP [orgn] AND alive";

   $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";

   echo $esearch.$SearchTerm;
}

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20 AND homo sapiens [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=HBB AND Escherichia coli [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=ABO AND Saccharomyces cerevisiae [orgn] AND alive


Is there a database of disease prevalence?

For rare diseases, you can use Orphanet. You can download an xml file with prevalence from orphadata.

There is an OWL version of Orphanet called ORDO, you can browse it on the EBI Ontology Lookup Service. See for example Arrhythmogenic right ventricular cardiomyopathy and point prevalence relationships in the lower right box. There are a variety of programmatic ways to access this, including the OLS web API, SPARQL queries, or downloading the OWL and using an ontology API.


What is the aim of insertion codes in the pdb file format?

When I wondered it myself years ago and I found a good explanation in a mailing list archive. I think it was this one, by Frances C. Bernstein:

  The problem with numbering started when people wanted to compare the
  'same' proteins from different species.  They found that there were
  the following possibilities that gave rise to differences:

  More or fewer residues at either end.
  Extra residues at various places within the chain.
  Fewer residues at various places within the chain.
  Different amino acids at the same place.

  Now imagine that residue PHE 195 is very important for the activity of
  the protein in species A.  But in species B it is residue PHE 197 and
  in species C it is PHE 212, because species B and C are not the same
  length as species A.

  Because people felt it was important to preserve the amino acid numbering for 'important' residues and to be able to readily discuss
  and compare the structures from different species, various people
  decided to try to number the proteins from species B, C, etc. to match
  the numbering used for species A.  In doing this, one must have gaps
  (missing numbers) where a sequence is shorter.  But what should one do
  when a sequence is longer?  This is the case where it is necessary to
  insert extra numbers and this is done by using insertion codes.

  Thus the insertion code is an integral part of the residue number and it is improper to ignore that field when using a PDB entry.  You
  must also allow for 'missing' numbers when using a PDB entry.

In the PDB format the insertion code is in column 27, following the sequence number (columns 23-26):

                      nnnnI
ATOM  11918  CZ  PHE D 100      -6.852  76.356 -23.289  1.00107.94           C
ATOM  11919  N   ARG D 100A     -9.676  74.726 -19.958  1.00105.71           N
...
ATOM  11970  CE  MET D 100H     -8.264  83.348 -19.494  1.00107.93           C
ATOM  11971  N   ASP D 101     -11.329  81.237 -14.804  1.00107.41           N

Together, columns 23-27 make a sequence ID, which is a mostly-numeric ID of the
residue position.

Note: It's not a full ID of the residue. When you have a point
mutation (a.k.a. microheterogeneity) you have 2+ residues with partial occupancy at the same position. To fully identify a residue both the sequence ID and residue name are needed.

As of 2020, the PDB (organization) uses PDBx/mmCIF as the primary format. The PDB format is only a legacy one. The presentation of the insertion codes (and many other things) in the mmCIF format are rather confusing, so it's worth explaining.

The mmCIF format introduced double identifiers for chains, residues and atoms: author's ID and "label" ID. Both look the same and are easy to confuse. For example: two chains named B and A by the author are renamed as, respectively, A and B by the wwPDB; both names are used in parallel: one chain is B(A), the other A(B). Poor design choice at every step. The same with the sequence ID. Originally, we had two tags:

_atom_site.auth_seq_id -- author's sequence ID, possibly with an insertion code (100A),
_atom_site.label_seq_id -- numeric-only sequence ID, numbered with respect to the sequence of this molecule (SEQRES).

Then, as I heard, numeric insertion codes were spotted in the wild. (In a database big enough any mistake that is not automatically checked for will happen). This made a problem with translation from the PDB format to mmCIF. The sequence number 100 with IC 1 would be interpreted as sequence number 1001. The numeric ICs must have been remediated in the meantime, because they absent as of 2020. But at that time (10-20 years ago?), instead of correcting the entries, the decision was made to move the insertion code to a separate column (pdbx_PDB_ins_code).
The documentation was never updated, so the auth_seq_id is still documented to store the sequence ID, not just the sequence number.

Here is an example from 3B9F:

                             !!  I                                      nn
ATOM   770  C CD2 . LEU B 2  46  ? -23.416 21.041  -27.401 1.00 19.94 ? 60   LEU H CD2 1
ATOM   771  N N   . TYR B 2  47  A -21.191 16.717  -28.119 1.00 18.73 ? 60   TYR H N   1

The sequence IDs here are 60 and 60A. Since the meaning of the insertion code is not widely known, it landed next to label_seq_id which has nothing to do with it. So again: the insertion code is just a part of the author's sequence ID (originally, in auth_seq_id), and the letter A above (10th field) is simply an extension to the number 60 (17th field).


How can I work out gram staining information from the species name?

The GOLD database can be used to retrieve a list of species names according to their Gram staining: go to the "Advanced search" page and select the either "gram-" or "gram+" under "Organism Fields" -> "Gram Stain". Then on the results page click on the number under "Organisms". On the results page you can select to have 1000 results per page, meaning that you would need to do some manual work.

I should add that as far as I know this information is based on what people submit when they register their genome projects, so I wouldn't expect this list to be curated or 100% accurate.

If you have access to the genome of the species you are interested in (you could download a representative strain for each species using ncbi-genome-download), you could use traitar, which can predict the gram staining, among other interesting broad phenotypes.


How can I use pandas agg here to avoid iteration?


circosJS should probably work well, based on your requirements: https://github.com/nicgirault/circosJS
I posted a short Observable notebook that shows what calls could be used to bring data from remote bigWig files into a circosJS plot:
https://observablehq.com/@alexpreynolds/rendering-bigwig-file-data-in-a-circular-representation
It is less intended as an explicit how-to in Javascript, than as a starting point for the end user with some of the function calls and data scale issues involved. Hopefully it can be useful to the reader for exploratory purposes.


Can I pilon-polish long reads with Illumina short reads to improve structural variant detection?

I'd recommend using Canu for correction rather than Pilon, because it has a component that is specifically designed for read-level correction. The newest version of Canu (v1.7) will track all reads through the correction and assembly process:

  Ensure that every raw read is either corrected or used as evidence for correcting some other raw read. This serves to rescue short plasmids in high coverage datasets, and it is no longer necessary to set corOutCoverage to achieve the same result.

More details here.


How to use variables while calling pybedtools

You have to take your variables and make them into a single string variable.

for example:

chr_ = "chr1"
start = 3000
end = 3402
as_str = ' '.join([chr_, str(start), str(end)])
# or
as_str2 = "{} {} {}".format(chr_, start, end)

Once you have the string built you can use the example in the bedTool docs.

a = BedTool(as_str, from_string=True)


Is there any difference between SNPs 'AG' and 'GA' in association analyses?

Could you please show us the context in which this appears, as you seem to be interpreting this differently to Devon.

If it's appearing as you say, GA and AG, then Devon is right, this usually means that the sequence goes ###AG### or ###GA###, which are two very different sequences.

If, however, as you're implying, the sequences are actually ###G### or ###A###, then semantically, the two ways of writing it have the same outcome. Normally, however, this is written as A/G or G/A.

Alternatively, you're looking at genotypes, which would usually appear as A|G and G|A. In this case, both individuals are heterozygotes, and in terms of the individual SNP they are the same. However the side of the pipe is relevant as this tells you which homologous chromosome it is on. This allows you to calculate LD with respect to neighbouring SNPs, and to determine the actual sequences of the proteins in that individual.


EMT/EMT-like processes in bioinformatics to study cancer progression

EMT describe cells' change in their state from being epithelial to the mesenchymal class. So, if a cancer cell line is gaining properties that allow it to move, it might reach the metastasis phase, were a cell can move to any point of the body and start a new tumor there.  As in the wikipedia page you linked:

  EMT has also been shown to occur in wound healing, in organ fibrosis and in the initiation of metastasis in cancer progression.

I am not aware of any kind of database for the EMT process, the most relevant resource might be the gene ontology, which has a term for it, which lists 185 (human) genes as related to this biological process.

ArrayExpress is a database of transcriptomics, it :

  stores data from high-throughput functional genomics experiments

You can evaluate the role of EMT in the data stored in the database.


How to merge sam files together with adding read groups

I have used https://github.com/ekg/bamaddrg to add read groups quickly to multiple sam files.  And then you can do a samtools merge of the tagged files.


How to I get the infix expression for the body in a SBML function defintion using python libsbml?

To do the translation you need to access the function definition. First, make sure there are function definitions in your model by using:

if model.getNumFunctionDefinitions() > 0:
   # continue

To get the infix of the body for the ith function definition use the following code:

f = model.getFunctionDeinfition(i)
body = f.getBody()
print (libsbml.formulaToL3String(body))

As a side note you can also get the jth argument for the function definition using:

f = model.getFunctionDeinfition(i)
arg = f.getArgument(j)
print (libsbml.formulaToL3String(arg))


Problems to Extract uniquely mapping reads from BWA MEM alignment

The issue is that you are filtering reads without regard to the pairing so you might filter out one of the mates if one aligns uniquely and the other does not.  You will need to filter out the un-mated reads from the file to have a valid BAM.

  As I do not see any XA or SA flags there anymore, I wonder what those can be?

I am not sure what this means.  The grep should remove those flags and only leave records without them.

As a side note if you don't want split alignments I would recommend using a different aligner such as BWA aln or BowTie2


making a customer reference by using cell ranger

Step 1:
Add the the GFP sequence to your fasta file

>GFP_name
ATG.....TAA

Step 2:
Annotate the gtf so that it contains all the necessary information to annotate your information

Add this line to your GTF file but modify where necessary:
"Make sure that all columns are tab separated"

GFP_name    AddedGene   exon    1   1966    .   +   1   gene_id "GFP"; transcript_id "GFP"; gene_name "GFP"; transcript_name "GFP"

The first column has to correspond to the header of the fasta sequence
start and end indicates the length of your sequence. I added a 1966 bp long sequence to my gtf

Code example:
step 1:

cat reference.fa gfp_seq.fa > new_reference.fa

step 2:

cat genes.gtf gfp.gtf > new_genes.gtf

-  gfp_seq.fa is your custom gene in fasta file
-  gfp.gtf is your custom gtf line or lines if you want to add more

You can also open your current reference fasta and gtf file in a text editor and copy paste the sequence and gtf line in there.

Then run cellranger mkref on your new_reference and new_genes file


Tree cut height for Consensus data input not working

Your data don't have obvious outliers so you don't have to worry about removing them. The problem is that both sets show prominent sample clusters. I would try to figure out what drives those clusters (is it technical or biological), and whether the variation should be removed. See the comments in WGCNA FAQ under point 5 (My data are heterogeneous. Can I still use WGCNA?). This looks like data from GEO, so look into the sample annotation for these data sets and see if any of the recorded sample characteristics align with the clusters you see in the sample trees.


Differential expression analysis for a subset of transcripts

It's safer to import everything.  You'll want all the data for normalization and dispersion estimates.


Selecting 65000 SNPs where AF is close to 0.5 in all or most populations

Just use bcftools view for filtering:

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf

To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.

$ bcftools view -h input.vcf.gz|wc -l
255

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf


Extracting specific sections from FASTQ determined by positions

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:

$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.

IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):

#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:

perl foo.pl file.fq 1,3,5-7


which NCBI tool is optimized to identify a species from a DNA fragment?

My conclusion from the tree is the fungi appears to be genuine and not an extant contaminant (which is very easily done). It is not Fomitopsis betulina because this is division Basidiomycota. The fungi here is division Ascomycota and these are very distinct from Basidiomycota.

There have been many examples of contamination in ancient DNA work, but this doesn't appear to be one of them. The rationale is given below.

The first fragment was not interesting and described 100% identity to ascomycetes. The tree above is a simple distance tree from the second fragment of 18S which shows it is an outgroup to the ascomycetes (green terminal taxa). The ancient fungi contains 3 (I think) unique nucleotide residues that have not been sequenced in any extant ascomycetes to date and sequences are unique to any known species. The location of the 3 unique nucleotides is tightly clustered within the locus with 100% identities either side. This is significant because it demonstrates it is not a PCR mosaic artefact (yep its been done in ancient DNA work and subsequently proven artifactual). In the latter PCR mosaic artefact I'm referring to (which I wouldn't cite), was a fungi/bacterial mosaic from memory.

Given the homology of 18S in the extant ascomycetes I would say this appears to be a new genus, however you would need to be sure this was not sequencing error and this is very old DNA, which is subject to degradation even when frozen. The implication of the tree is quite unusual, it implies widespread extinction of this group of fungi in recent history (5000 years is slight in evolutionary terms) BUT it could be a sampling bias, i.e. no one has sequenced those fungi and NCBI is not representative.

I assume that you have performed phenotypic/morphological analysis on the spore to believe it is Basidiomycota.

The fragment in yellow is the query and the sequence next to it was the sequence already deposited on Genbank of this ancient isolate. You would need to look up the biology and ecology of ascomycete to understand its significance, but it is in a lineage of its own, so maybe difficult to interpret by comparison to extant species. You could concatenate the two alignments based on each fragment, but the problem with 18S is that is a hard to model and you'd need to perform secondary RNA structure modelling to assess whether these mutations would be viable.

In conclusion, it definitely looks genuine but it is not the species you think it is, in fact not even close.


Can I ask DESeq2 which variable alone explains which gene's behavior the best?

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


Order of batch effects removal, data imputation and library size normalization in scRNA-seq data

MAGIC assumes input data has been both library-size normalized, and either log or sqrt transformed prior to imputation (see also: MAGIC tutorial). Additionally, any graph-based methods (MAGIC, PHATE, t-SNE, UMAP, spectral clustering, Louvain, etc etc) will give flawed results if your data contains a batch effect, since the neighbourhood graph would reflect that structure of your batch effect, and worse, imputation would further reinforce this batch effect.

Thus I would recommend the following pipeline:

Library-size normalization
Square root (or log) transform
Batch effect removal
Imputation

Regarding your update, the reason you don't see the batch effect in the raw counts data is simply that the batch effect is not visible in the most highly expressed genes. Prior to transformation, the principal source of variation in your data is simply the expression of the most highly expressed genes - this is essentially masking the source of the batch effect, not removing it. I recommend never working with raw molecule counts in scRNAseq, as the raw counts data hides much of the heterogenity in your dataset, which is precisely what you are looking in doing single-cell RNA-seq.


How can I translate infix to mathml using the python interface to libsbml?

It’s quite straightforward to convert infix to mathml using libsbml. Make sure you have libsbml installed in python by typing

pip install python-libsbml

at the windows/Mac/Linux command line or

!pip install python-libsbml

if you are running a Jupyter notebook. Some tools such as Tellurium already have libsbml preinstalled.

At the python command line type

import libsbml

to import the package. To enter the infix expression type the following:

p = libsbml.parseFormula ('k1*s1')

To generate the mathml type

print (libsbml.writeMathMLToString (p))


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).


Removing repeated reads from nanopore 1D² reads

You can exclude the 1D² reads by generate a list of the ids from the 1D² reads, then filtering those reads out from the 1D fastq files. Here's one way to do that using some accessory scripts that I have written:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | awk '{print $2}' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

In order to deal with the concatenated read IDs, it's necessary to switch to a language like Perl that supports the use of backreferences for the splitting. The following perl command (replacing the first line in the above procedure) inserts a line break after the 36th character in the 2nd value, then prints that value:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt

So a full sequence of operations would look something like this:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq
$ cat 1d2_reads.fq filtered_1d_reads.fq > unique_reads.fq