Transform a R data.frame into a krona plot without krona tools?

First download any Krona plot (.html) file or make one by using the Excel sheet from the Krona Package. Open in text-editor and remove anything below "datasets" tag and/or the first "node" tag till the end of the document.

Here is the path to go: 1) get your data.frame in hierachical order 2) transform to data.tree object 3) transform to xml file 4) gsub until all the tags match krona annotation

options(scipen=999)
options(stringsAsFactors=FALSE)
library(data.tree)
library(xml2)
library(magrittr)

#get your data into a hieracial format in a data.frame
t <- structure(list(freq = c(34L, 127L, 42L, 123L, 343L, 13L, 104L,
21L, 135L, 91L, 167L, 148L, 90L, 46L, 25L, 79L, 1726L), king.n = c("",
"Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Plantae",
"Plantae", "Plantae", "Plantae", "Plantae", "Plantae", "Plantae",
"Plantae", "Plantae", "Plantae"), class.n = c("", "Dothideomycetes",
"Dothideomycetes", "Dothideomycetes", "Dothideomycetes",     "Dothideomycetes",
"Dothideomycetes", "Liliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida"), fam.n = c("",
"", "Cladosporiaceae", "Didymellaceae", "Pleosporaceae", "Pleosporaceae",
"Pleosporaceae", "Poaceae", "", "Asteraceae", "Asteraceae", "Asteraceae",
"Asteraceae", "Fabaceae", "Hydrophyllaceae", "Hypericaceae",
"Malvaceae"), gen.n = c("", "", "Cladosporium", "", "Alternaria",
"Pyrenophora", "Pyrenophora", "Phleum", "", "", "Achillea", "Centaurea",
"Jurinea", "Trifolium", "Phacelia", "Hypericum", "Malva"), spec.n = c("",
"", "cladosporioides", "", "", "", "tritici-repentis", "pratense",
"", "", "", "cyanus", "", "", "", "", "")), .Names = c("freq",
"king.n", "class.n", "fam.n", "gen.n", "spec.n"), row.names = c(1L,
3L, 4L, 5L, 8L, 9L, 11L, 15L, 16L, 17L, 18L, 19L, 21L, 22L, 25L,
26L, 28L), class = "data.frame")

#if any sequences can not be assigned to the highest hieracial level give them a random name like "NA"
t[t[,"king.n"]=="","king.n"] <- "NA"
#define the hieracy (do not include the counting var here), the first space is the root category
t$pathString <- paste(" ",t$king.n,t$class.n,t$fam.n,t$gen.n,t$spec.n,sep = "/")
#transform to a data.tree
t <- as.Node(t)
#define the sequence frequency calculations
t$Do(function(x) x$freq <- Aggregate(x, "freq", sum), traversal = "post-order")
print(t,"freq")
#get all frequencies for each node
repval <- t$Get("freq")

#tranform into a nested list and afterwards into a xml document to get the nesting
t <- t %>% as.list(.,mode="simple",unname=T)
t <- as_xml_document(list(root=t))%>%as.character

#from here it is just a bunch of regex
#split the xml file into one line for each <
temp <- strsplit(t,"\n")[[1]][2] %>% strsplit(.,"(?=<)",perl=T) %>% .[[1]]
temp <- paste0(temp[c(T,F)],temp[c(F,T)]) %>% gsub("[0-9]","",.)
#remove the names nodes and replace by krona annotation
#replace the closing named tags by generic closing "node" tags
temp <- gsub("</(.*)?>","</node>",temp) %>% gsub(">.*$",">",.)
#replace the named opening tags by generic "node" tags with the "name" attribute, replace all (because they are incomplete) numeric values with a placeholder
temp <- gsub("<(?=[^/])(.*)?>",
     '<node name="\\1"><magnitude><val>__\\1</val></magnitude><score><val>1</val></score>',temp,perl=T)
temp %<>% gsub("root"," ",.)
#find the placeholders to insert the frequencies and execute
repos <- grep("__",temp)
for (i in repos)
  temp[i] <- gsub("^(.*)(__.*?)(?=<)(.*)",paste0("\\1",repval[which(repos==i)],"\\3"),temp[i],perl=T)

#add the end of the .html file
temp <- c(temp,"</krona>
","","")  %>% unlist
tout <- paste(temp,collapse="\n")

#append your data to the Krona plot
write(tout,file="Krona1.html",append=TRUE)

Et voilá!


How to quantify similarity of genomes and find differences in set of S aureus genomes?

... stuff of original post deleted.
On second thoughts what you might be doing is a template assembly of your genomes. It is a possible interpretation of the 1.2. fasta sequences above (i.e. 1 is template).

A microbial genomics professional would advise also performing a de novo assembly, particularly if you are interested in presence and absence. The reason is that if a gene is present in your query and absent in your template, you will miss it. Again it all depends on what bacteria are being assessed, some are more prone to "genetic islands" than others.

You need a collaborator beyond that.

For the onlookers here bacterial genetic behaviour is very different to eukaryotes and what they get up, switching DNA etc.., to would appear bizzare from a higher eukaryotic world.

You mentioned a comparative analysis of 500 strains: I have worked on this in MRSA. Note, the bacteria is important considering the approach you will adopt.

Anyway, you want a single aligned file to produce a phylogeny, Bayes or likelihood. This is a model of point mutations. Bootstrapping doesn't really help because of SNP differences between isolates.

There is a complex problem. Generic phylogeny can hit the buffers in my opinion because often:

low numbers of SNP differences between isolates.
The error across the genome is poorly defined
Multi-clonal infection ignored and your isolates not cloned.
The other issue is the QC of the genome, I've been fairly astonished by the variation (it affects the tree)

Put it all together and you can get a nice tree but the topology of a given MLST is in my opinion unlikely to be correct. A referee may not bother with this, but they might.

You then map the phenotype against the phylogeny (the tips of the tree) and look for clusters. HOWEVER, there are problems:

Drug-resistant e.g. mec cassette (methicillin-resistant ) are often on plasmids and these are lost during isolation. Albeit they can integrate onto the genome. So an isolate can be drug resistance but you fail to find the gene.
Drug-resistance transmission is unlikely to cluster, its too quick, so you don't see nice tight clusters against a theoretically perfect tree
The best approach is drug-resistance in culture.

Hierarchical clustering, has many meanings. It is used in presence/absence data, I demonstrated the method failed (on bacteria) based on cluster analysis. It is currently revised and used in unsupervised deep (or machine) learning prior a training method.  The question then is what is what are you modelling? Presence/absence (of genes), point mutations, epidemiological data, drug-resistance/non-drug resistance?

My assessment is that given the potential complexity of these biological scenarios and the quality of the question you should seek formal collaboration and this is before considering what bacteria are being assessed. Some of epidemic drug-resistent bacteria rip up this rule book.


the model matrix is not full rank : This is a classic question which a biologist face without clear understanding of the model design

> table(doldat$Group, coldat$Mutation)

        DNMT3A-R882H FLT3-ITD IDH2-R140Q KRAS-G13D MAU2-Q98H No NPM1-L287ins(TCTG) NRAS-G13D STAG2-R614 TET2-E1357 TET2-R550
  HSC              0        0          0         0         0  7                  0         0          0          0         0
  Blast            2        4          2         1         1  0                  4         1          1          6         1
  CMP              0        0          0         0         0  8                  0         0          0          0         0
  GMP              0        0          0         0         0  7                  0         0          0          0         0
  LSC              2        3          0         0         1  0                  0         1          1          0         0
  Mono             0        0          0         0         0  6                  0         0          0          0         0

Look at this table, you see by eye how some mutations are nested with Group. If you want to include it then you will need to subset the experiment into multiple separate runs, depending on the question you want to answer.


So B-factors, for "blur factors", are how much movement an atom has in the crystal structure. They are not reliable. B-factor putty in PyMOL is the most you should push them.
B-factors bad
This is partly a result of flexibility, but it is highly influenced by crystal packing, which is not natural. Many enzymes, for example, are seen with an active site entrance way that is closed up in a crystal structure whereas in an MD run there is a gate loop that flops around.
Additionally, a B-factor is a single value that masks lots of forms of movement. Partial occupancy after all is at the discretion of the experimenter, so you might even have a residue in two or three low energy states that result nevertheless in B-factor disorder.
Furthermore, incorrectly fitted residues will have a higher B-factor.
Also, the meaning of B-factor are mighty cryptic, and depend on the refinement —the units are Å^2 as opposed to say Å, like a standard deviation.
MD RMSF good
The best way to determine the movement of a residue is to do an MD run. Gromacs has very nice tutorial and is the most accessible.
The metric in that case is the RMSF, which is the RMSD relative to the average. Gromacs has a utility for it called rmsf. Alternatively RMSD relative to your original structure is fine —but less of a metric of movement.
MD is complained about for three reasons, none are valid here:

Tricky to set up: this is a straightforward case
Lengthy: You would not be after large scale conformational changes so a limited amount of steps would suffice —possibly even on a laptop on a single core for a couple of hours (cluster is always better though!). Running at different temperatures and switching the intermediate steps (replica exchange) is often done to speed things up, but would just confuse the results.
Small molecules: A nanobody has no cofactors or small molecule ligands so you'd avoid the topology parameterisation issues. If Eukaryotic made, your nanobodies may be glycosylated — play it safe and ignore.

B-factor extraction
If you want to extract the B-factors, you can use pymol as a python module —installed with conda and not the binary/exe installer. And use python as follows:
pdb_filename = ""
selector = "name CA"
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(pdb_filename)
     bfactors = {atom.resi: atom.b for atom in pymol.cmd.get_model(selector).atom}

From there do whatever you fancy with the b-factors. For alternative selector strings see Pymol selection algebra. One operation is to divide the values by 8 pi^2 and take the square root of that, which is the amplitude of oscillation (U) of the atom around its equilibrium position —ie. that stated in the coordinates.


DNA bioinformatic search and methylation analysis by qMSP

It's likely that they accomplished this by using the UCSC Table Browser to:

Find the loci of the two isoforms you mention
Finding CpG islands (Group = Regulation, Track = CpG Islands) near those transcription start sites
Design primers using Applied Biosystem's software

Alternatively, they list the exact primers used in Table 1.

They don't detail exactly how they did it, so if you have specific questions it might be worth it to email the corresponding author


The data are imputed using values from the same column

Here's a simplified example to show what's going on:

A <-c(1,2,3,4,5)
B <-c(2,3,NA, 2,4)
C <- c(2,3,4, 6, 1)
data<-data.frame(A,B,C)
data

#>  A  B C
#>1 1  2 2
#>2 2  3 3
#>3 3 NA 4
#>4 4  2 6
#>5 5  4 1

# I'm using scale = F to make things simpler
knnImputation(data, k = 3, scale = F, meth = "weighAvg", distData = NULL)

#>  A        B C
#>1 1 2.000000 2
#>2 2 3.000000 3
#>3 3 2.594272 4
#>4 4 2.000000 6
#>5 5 4.000000 1

# calculate euclidean distances, dropping column of interest
dist(data[,-2])

#>1        2        3        4
#>2 1.414214
#>3 2.828427 1.414214
#>4 5.000000 3.605551 2.236068
#>5 4.123106 3.605551 3.605551 5.099020

# use k rows with smallest distance to row of interest to calculate weighted average
# using exp(-dist) as weights
weight_sum <- weight_sum = exp(-2.828427) + exp(-1.414214) + exp(-2.236068)
2*exp(-2.828427)/weight_sum + 3*exp(-1.414214)/weight_sum + 2*exp(-2.236068)/weight_sum

#>2.594271

If you want to impute using values from the same row, all you have to do is transpose the matrix:

t(knnImputation(t(data), k = 2, scale = F, meth = "weighAvg", distData = NULL))

#>     A        B C
#>[1,] 1 2.000000 2
#>[2,] 2 3.000000 3
#>[3,] 3 3.086729 4
#>[4,] 4 2.000000 6
#>[5,] 5 4.000000 1

(I had to lower k because there are only 2 extra columns)


Your question is not clear, you seem to mix the relevant terminology...

First, you cannot run WGCNA for individual samples, only for individual conditions (with or without WT, depending on what question you're trying to answer and assuming you have enough samples in each condition). If that's what you have done, fine. The number of modules in each independent analysis will generally be different, that's because each condition is different. The module identification built into WGCNA as the default approach does not allow the user to specify the number of modules beforehand, but if you insist on getting a pre-specified number of modules, you could run WGCNA using the "manual" way (for example, Tutorial I, section 2b) and substitute another clustering method (e.g., PAM) that allows you to specify the number of clusters. You may find the results inferior though.

Second, module preservation does not result in new modules. Module preservation evaluates whether modules found in data set A are preserved in a different data set B. To get new, common modules, you could run a consensus network analysis on pairs of conditions or all 3 conditions together. The number of modules you would find in a consensus analysis cannot be predicted from the number of modules you find in individual data sets.


The dbNSFP plugin from VEP accesses tables of data for each variant from dbNSFP and pulls out the values. dbNSFP provide their SIFT scores in that format: a score for every transcript affected by the variant, all on one line. The lookup is just for the variant, not for the variant/transcript combo, so they provide scores for all variant/transcript combos. You can also get a column that gives you a list of the transcripts or proteins (Ensembl_transcriptid or Ensembl_proteinid) in order so you know which score goes with which transcript.

A better way to get SIFT scores with VEP is to get them directly from VEP, rather than using dbNSFP. This will get the SIFT score that goes with the transcript on the line with the relevant transcript.


Difference between clinical significance and statistical significance

This is non-parametric statistics, the mean requires confirmation to the normal distribution (but there are exceptions), or perhaps better put it requires a fixed relationship between mean and standard deviation. Non-parametric stats always use the median and are distribution free. Wilcoxon test/ Mann-Whitney U test  are used to replace t-test. The Kruscal-Wallace test replaces the 1-way ANOVA
Finally the Friedman test replaces 2-way ANOVA but has its critiques.  The t- and ANOVA tests require normality, F-max is one approach to determine normality.


You need to run Bracken to the Kraken2 report output to estimate abundance.


Can Maven (mass spectrometry software) read mzML files?

No, Maven still can read only mzXML files.


Generating burrows-wheeler-transform and a suffix array of a DNA sequence with less memory

No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the files.


Program Quality control for data sequencing (nanopore)

A simple start would be to use MinIONQC, which is also available via bioconda under the name r-minionqc.


Rdkit.Chem.Draw.savefig() creates empty figures

For some reason you have to set bbox_inches to tight. Without that rdkit does not fit the molecule into the canvas and might plot an empty figure.

fig.savefig('mol.jpeg', bbox_inches='tight')


Usually it doesn't actually matter, but if you want it to be very correct:

RGID=SomeLibraryID \
RGLB=SomeLibraryID \
RGPL=illumina \
RGPU=Some_Sequencer_ID \
RGSM=SomeSampleID

That assumes you have one library per sample and aren't splitting the output BAM files by lane.


Col-0 is, to my knowledge, the standard reference Arabidopsis thaliana ecotype, so the easiest way to get the sequences for all genes is to just download them from Ensembl. The nucleotide sequences are available here and the protein sequences here.


Get & Annotation of isomiR (isoform of miRNA) profile from RNA-seq

The obvious answer is to break up the file size and upload each separately,
First I'd check the file size,

wc infile.fa

Output

  XXXXX      YYYY     ZZZZZ infile.fa

The number you want is XXXX - the number of lines in the file.

You'd then decide how many files you wanted to break it up into ... in your case divide XXXXX by 3. Then .. for example

sed -n '1,5000000p' infile.fa > outfile.fa
sed -n '5000001,2000000p' infile.fa > outfile2.fa
sed -n '2000001,4000000p' infile.fa > outfile3.fa

I'm not certain this would be cool however. An miRNA expert could answer this question very quickly. I think is you should check your data for sequence duplication. Simply piping the fasta file through a simple distance matrix and screening for "0" values would confirm this.

Creating a frequency table of duplicate RNA must be a feature of Bioconductor, but in any case screening for duplication (technically duplicate records) is a formal SQL database command.

It can be done a Perl/Python but its a bit clunky for a data set this size. I'd do it by setting the sequence as the key within a hash, the inefficiency has been solved (but its a weird way to do it).

You would have a database for "headers" for a given sequence, which would take up space, but I would assess the potential level of duplication first.


Automating download of results of Drug and Drug Target Mapping

Simpler option would be to not use the query interface  and directly parse the CSV file found on that page.


As Devon said, removing the which will give you W as a logical vector that can be used to index the sorted p-values. To get a logical vector that can be used to subset the original (unsorted) values you’ll need to go a step further:

candidates = is.element(1 : L, order(adjusted.p.values)[W])

However, I’m not sure this will work with your data: the name adjusted.p.values suggests that these values are already FDR-adjusted. So you cannot/should not perform additional correction on them. If they are already adjusted, then you can get candidate loci simply by

candidates = adjusted.p.values < q


How to input the opposite strand of DNA?

Yes, you're going to want to take the reverse complement of your string in order to analyze the other strand, that is the standard.
Also, it is worth pointing out that you may want to use a generator expression instead of a list expression for your chunk_string function:
def chunk_string(s, n):
    return (s[i:i+n] for i in range(len(s)-n+1))

This will result in a generator that yields kmers as the Counter constructor iterates through them, which means that you don't need to first construct and store the entire list of kmers first and then pass to Counter but instead can feed them directly to it. Not essential for small strings, of course, but a good practice and useful for larger strings. See this SO post for a better explanation.
You can also use Biopython's reverse_complement function in order to easily obtain the reverse complement.


What does the number mean in an HGVSp annotation?

It's the position in the protein. p.Arg452Pro means that there is a variant that changes an Arginine residue at position 452 of the protein into a Proline.

And yes, it is absolutely possible to have the same HGVS p. notation on the same gene but on different transcripts. More than possible, it is quite common. Your particular example isn't actually a real variant since that position on chromosome 4 is intronic in both hg19 and hg38, but TP53:p.R175L1 is a real example of what you describe:

As you can see in the image or by visiting the link, this variant on TP53 has the same HGVS p. notation on 6 transcripts. It also has a different p. notation on two more transcripts that are not shown in the image as well as no p. notation for a few other transcripts in which it doesn't affect the protein product.

1Disclaimer: The link goes to https://varsome.com, a free human variant search engine that is developed by the company I work for.


Choosing Fisher's Exact or Binomial test for overrepresentation in PANTHER

For a guide on using the tests in PANTHER, see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6519457/.
In answer to this specific question, the defaults are the best choice (Fisher's Exact Test, with FDR correction). Fisher's Exact test gives more accurate statistics for relatively small counts (e.g. when only a few genes in your list are associated with a particular GO term). The FDR is a better correction to use for these types of analysis in general, especially given the large overlap between classes.
The other options are there mostly to support replication of older analyses, e.g. to replicate a published analysis that used the PANTHER over-representation tool prior to 2018, when the binomial with Bonferroni correction was the default.


Fast filtering of intervals not falling within a certain distance from known genes

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

Convert GTF annotations to a sorted BED file of genes:

$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

Sort your intervals, if unsorted:

$ sort-bed intervals.unsorted.bed > intervals.bed

Filter intervals against 1kb-padded genes:

$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed

We first --merge genes to simplify dealing with overlapping genes and further speed downstream filtering.

The file intervals.filtered.bed will contain intervals from intervals.bed which lie outside of a 1kb-window around merged gene space.

To run all of this as a one-liner (assuming use of the bash shell):

$ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed

This creates the same output result, but avoids creation of intermediate files and so should run even faster.

Note: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using --range takes care of ranges that would have a negative start position, setting their start coordinates to 0. Bounds aren't really relevant here, in any case, because the bedops --not-element-of step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.


In Uniprot how do I know which variants are mutations and which are polymorphic?

This answer is very late but may still help others:
Have a look at the UniProt index of human polymorphisms and disease mutations at https://www.uniprot.org/docs/humsavar

Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.


Something like vcf-isec? Quoting:

  Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.


Way to get genomic sequences at given coordinates without downloading fasta files of whole chromosomes/genomes first?

using a http request.

if there is a DAS server, you can always use this protocol to download the xml -> fasta. see https://www.biostars.org/p/56/

$ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200"  | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n'
ggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctacta

or use the UCSC utility twoBitToFa which works with remote files.

If you have an NCBI accession number you can use the E-utilities with seq_start and seq_end.

$ wget -q  -O -  "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AE014134.1&seq_start=100&seq_stop=300&rettype=fasta"
>AE014134.1:100-300 Drosophila melanogaster chromosome 2L complete sequence
TGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGATGATAATATA
TTCAAGTTGCCGCTAATCAGAAATAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATG
CGAGAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG


bowtie2 options when mapping stranded single end reads

I repeat what user172818 says...don't waste your time on bowtie or tophat.  Use STAR or HISAT2, or a pseudoaligner like kallisto with a transcriptome.


Finding the members in a confusion matrix

This does not look a great confusion matrix, but it is a very, very cool approach. The cool thing about machine learning is the ability to examine lots of models within confusion matrices.

To answer your question, you want a heat map of the regression weights to establish who is who, i.e. which patient is false positive etc... This is the only way I'm aware of to analyse a neural network. This looks like R, I only know Scikit learn.


Database for proteome-wide predictions of protein structures

There are three initiatives I know of to have a go at this:

PConsFam, which collects data from this paper.
The Baker group's metagenomic study which you mention in the question.
The recent DMPfold work from our lab. This compares its results to the above two studies and discusses the effect on various model organisms.

The second and third of these don't provide a database per se, but the models are available for download and are linked to the Pfam database.


How to efficiently get human gene names from NCBI based on a large list of SNPs

Modify SNPs and gene annotations for your genome assembly of interest. This example shows how one might do this for hg19, for instance:

Download SNPs and convert to BED:

$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz \
    | gunzip -c \
    | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - \
    > hg19.snp151.bed

Filter this BED file for entries of interest:

$ grep -wFf snps-of-interest.txt hg19.snp151.bed | cut -f1-6 > hg19.snp151.filtered.bed6

Download gene annotations and convert to BED:

$ wget -O - ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
    | gunzip -c \
    | awk '($3 == "gene")' \
    | gtf2bed \
    > gencode.v19.genes.bed

Note: In the future, newer versions of Gencode annotations may be available for your reference assembly. Always check their site to confirm what is most recent.
Map genes to SNPs:

$ bedmap --echo --echo-map-id hg19.snp151.filtered.bed6 gencode.v19.genes.bed > snps_with_associated_gene_names.bed

A list of genes that overlap a SNP of interest will be in the seventh column of the output.

If this must be done within Python, you could use something like subprocess.call(...), where ... includes commands wrapped in bash -c.

I don't recommend this because of the need to escape quote marks and deal with Python's quirks in running command-line tasks, which makes a Python-based approach very difficult to set up and maintain. I'd suggest learning some basic bash scripting to make this easy and fast.


If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?

I actually disagree with that.. I guess I write this down as a discussion.
The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.
However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:
plot(dpois(0:10,lambda=1),ylab="Probability",xlab="No of reads")

The probabilities are like:
data.frame(n=0:5,p=dpois(0:5,lambda=1))

  n           p
1 0 0.367879441
2 1 0.367879441
3 2 0.183939721
4 3 0.061313240
5 4 0.015328310
6 5 0.003065662

The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206


When working with Bioconductor and the curatedOvarianData dataset, what are the unit values of the expression data?

These values represent the expression of a particular gene in each sample. While expression of an immune system gene like CXLC12 might be associated with the aggressiveness of a cancer, that's not going to be the case for the majority of genes.

Looking at the numbers, these look like they are probably log$_2$-transformed values. In order to know for sure, plotting a density distribution of the dataset is a good idea:

plot(density(exprs(TCGA_eset)));

If the values are in the range of around 0..15, then it'll be log$_2$-transformed values, whereas if they're in the range of 0..100000, it'll be "raw" expression values. The log-transformed values are preferred for data analysis because gene expression tends to be close to a normal distribution in log-space.

These numbers are best treated as relative counts (i.e. the most appropriate unit is "expression"), as variation in the preparation of the sample (or the sequencer, or the microarray) can influence the results that come out the other end. Converting them to something more meaningful (e.g. "transcripts per cell") requires normalising against a known reference sample.


What pitfalls exist with running FastQC on a bam file?

I don't see why you are trying to compare two totally different measurements.
When FASTQC says a read is unique, it means that no other read shares its sequence.
When STAR says a read is uniquely mapped, that means it maps to one and only place.  That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.
You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.
FASTQC was designed for DNA sequencing, not RNASeq.  Some of its warnings make no sense in the context of RNASeq.  If I was doing RNASeq of a filtered set of pancreatic alpha cells,  I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.
You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes.  The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.
That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong.  Find out where all those reads are aligning, you might have rRNA contamination.  If you are doing an organism with a high level of duplication in its genome that might be normal.


Rearrange your text file into a sorted BED file:

$ tail -n+2 segmentation.txt | awk -vFS="\t" -vOFS="\t" '{ print $2, $3, $4, $1 }' | sort-bed - > segmentation.bed

Use BEDOPS bedmap with --echo-map-id-uniq and --count operators on merged regions:

$ bedmap --delim '\t' --echo --count --echo-map-id-uniq --fraction-map 1 <(bedops --merge segmentation.bed) segmentation.bed > answer.bed

The file answer.bed will contain merged intervals from segmentation.bed in the first three columns, the number of intervals from segmentation.bed which are contained entirely within the merged interval in the fourth column, and a unique listing of IDs of those contained intervals in the fifth column.


Multiway intersection of 8 lists of genes

EDIT: Sorry, I did not see the R tag and wrote the code in bash. I still left the bash example below and added an R version.

Depends on what you want to do, but there are two common ways in which you could store this. Either as wide format:

      h2  h4  h6 ...
Gene1 1   1   0
Gene2 0   1   1
Gene3 1   1   1
...

Where 1 denotes presence and 0 absence, or as long format:

      timepoint
Gene1 h2
Gene1 h4
Gene2 h4
Gene2 h6
Gene3 h2
Gene3 h4
Gene3 h6
...

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

Generate table in wide format:

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

R:

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

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

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

Bash:

Note: you will might need to add headers to this table

for timefile in $yourfiles;do
  # make a file with all gene names
  awk -v tp=$timefile '{print $1}' >> names_file
done

sort names_file | uniq > tmp && mv tmp names_file

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

Generate table in long format:

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

R:

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

Bash:

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

EDIT2: To follow up on this, here is a simple way (using the package dplyr) to count the number of common genes between all possible combination of time points. Assuming you have your data in long format:

df <- data.frame(genes=c("G1","G1","G2","G3","G3"), tp=c("h2","h4","h2","h4","h6"))

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

You can get the co-occurrences using :

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

And visualise it using for example a heatmap:
(ofc it looks stupid because I have a tiny sample where not all combinations are even present)

# Visualize number of genes present in each combination
library(ggplot2)
ggplot(df, aes(x=tp.x, y=tp.y, fill=N_genes)) + geom_tile()


How to plot p-values in a circular barplot?

The adjustcolor function from grDevices can be used to generate colours with transparency. What you want is colours based on the pathway and transparency to form a colour gradient for the p-value.

adjustcolor( "red", alpha.f = 0.2)

For the ggplot example, the aes function has an alpha parameter.

ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found), alpha = Entities_pValue) + ######################################
geom_bar(aes(fill = Entities_found, alpha = Entities_pValue), stat = "identity", position = "stack") +
coord_polar() +
theme_minimal() +
 guides(fill = FALSE)


about the cpm normalization after using normalization factor

In fact the cpm function will use the size factors (TMM) that were calculated with calcNormFactors. In the absence of the size factors cpm would do a naive per-million scaling so only correct for differences in library size. The power of the size factors is to also correct for library composition. I recommend this video for details. In short: What you do is fine. I would suggest though that you use the log = TRUE option for cpm as one typically performs downstream applications (even when using the Z-score) on the log scale. Or alternatively t(scale(t(log2(z+1)))).


Building network for biomarker interaction analysis

You can start by looking at the correlation (if data is non-parametric) of each variable with each variable of each data point you have. One thing that might change is the correlation between variables along the disease.

Next you can determine which relationships are constant along the disease progression/clinical trial. Identify also the ones that don't show a constant relationship. To do so you can use WGCNA (although with just 10 variables is too low to expect a scale free network you can use to correlate the variables and group them by correlation profile)

Then the question is which variables are more important, those who kept the same relationship or those which change?
Use the eigengenes of those variables (separately by those that have a constant relationship and those which don't) to estimate the dependency with the parametric variables. Or you simply can do an ANOVA of all the variables with each non-parametric variable.


What is the difference between a MeSH concept and a MeSH descriptor

MeSH is organized into Descriptor > Concept > Term hierarchy, where all the synonymous <term>'s are grouped together under a <concept>. For example, consider the MeSH terms,

Blood-Brain Barrier
Hemato-Encephalic Barrier
Brain-Blood Barrier

All these terms are group under the umbrella concept of "Blood-brain barrier".


What kind of "gff" format does bioawk parse?

I would consider the description there a bug. The filter is actually the strand, strand is the frame, group is the attribute, and attribute does nothing. These are really meant to be the 9 columns.

Edit: There's a bug report related to this.

Edit 2: I've made a pull request to clarify this and fix the aforementioned bug report.

Edit 3: I realized that I never directly answered the title of your question (mea culpa). bioawk itself will work with gff, gff3, or gtf files. It really is just treating them as tab-separated files with named columns (this is surprisingly convenient, since it's a PITA to remember what column does what).

Edit 4: The PR has been merged. If you install from github then you'll see corrected field names.


Greedy Motif Search Using Probability Matrices

Assuming you're talking about Implement GreedyMotifSearch from Rosalind, while I haven't actually finished this particular one, I've got an idea why TTC,ATC,TTC,ATC,TTC doesn't work. Here's my explanation from the "Questions" section of that problem:

  The profile is created from the first kmers in each string, and the most likely string from those is CAG. The answers will generally be similar to this.

I've been trying to do Rosalind by only using the information provided in the questions; this is one of the many problems where there wasn't enough information in the Question for me to arrive at a correct solution when I tried it. It's possible that the debug datasets were added to this question since I last looked at it, as they seem to describe this situation:

TEST DATASET 1:

Input:

3 4
GCCCAA
GGCCTG
AACCTA
TTCCTT

Output:

GCC
GCC
AAC
TTC

  This dataset checks that your code always picks the first-occurring Profile-most Probable k-mer in a given sequence of Dna. In the first sequence (“GCCCAA”), “GCC” and “CCA” are both Profile-most Probable k-mers. However, you must return “GCC” since it occurs earlier than “CCA”. Thus, if the first sequence of your output is “CCA”, this test case fails your code.


Finding discriminating positions for a sequence alignment column

Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.


How to get all PDB homologs from Uniprot (mapping + BLAST)?

For the correspondence between PDB and Uniprot entries you can use SIFTS -- a semi-automated mapping between PDB and UniProt maintained by PDBe.

The pipeline that creates the mappings uses BLAST and a few other criteria to decide which UniProt entry should be assigned to each PDB entry. The SIFTS website has all the data in CSV files.


This doesn't use R, but a simple and quite convenient method is to use bedtools getfasta with your fasta file. You'll need to make a BED file out of the RepeatMasker output:

awk 'BEGIN{OFS="\t"}{name=sprintf("%s_%s", $10, $11); strand="+"; if($9 == "C"){strand="-"}; print $5,$6-1,$7,name,".", strand}' repeatmasker.output > rmsk.bed

Then, given that and your genome.fa fasta:

bedtools getfasta -name -s -fi genome.fa -bed rmsk.bed -fo output.fa

This will produce a multiline fasta file names like AluSx_SINE/Alu (due to the -name option) where the sequences are reverse complemented if appropriate.

You can install bedtools either via conda (conda install -c bioconda bedtools) or github.


Autodetect max number of cores and pass as an argument in Nextflow

I have not read the full Nextflow manual in depth, but I think there is no Nextflow command to detect the maximum number of threads you can use.

However, you can try to automatize this task by using the nproc command.

nproc gives you the number of processing units available. An idea would be to include something such nproc - 1 as the max number of threads for Roary.


If you want to download a file, you need to use the link to that file. Your original attempt, wget  https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474 wouldn't work since that's a link to the trace page of the relevant run. If you want to download something else, just right click on the link (the one in the screenshot in your question), copy the URL and use that:

wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam


Cuffmerge: EOF marker is absent. The input is probably truncated

There is a post about this on seqanswers. The poster says the following about cuffmerge behavior:

  I kept the temporary files while running cuffmerge. It builgs a SAM file from transcripts.gtf files. Then it seems first it reads it as BAM since I read samtools take the files as BAM default. Then it finds it is not a BAM so it tries SAM. Then it can read it and resume the task.

So it sounds like the BAM file is generated from GTF as a kludge inside cuffmerge, not the other way around.

Did you change any of the intermediate files or delete anything? Is there maybe a different cuffmerge run that it might be colliding with?

Did you get any output? it sounds like this error is frequently observed when actually the command completes fine, e.g. https://www.biostars.org/p/189139/. See other posts by googling the problem maybe, e.g. https://www.biostars.org/p/344439/

(I know this isn't a proper answer but it was too involved for a comment.)


bedtools: get name column of alignment file A

I would recommend using bedtools intersect for that type of comparison.  This will allow you to see the specific segments that overlap.  bedtools coverage is used more for determining how many things cover a region but not which ones.

The documentation linked above will provide a much better desciption of intersect than I would be able to in this post but I would recommend using the -loj or -u depending on what you want.

-loj    Perform a “left outer join”. That is, for each feature in A report each overlap with B. If no overlaps are found, report a NULL feature for B.

-u  Write original A entry once if any overlaps found in B. In other words, just report the fact at least one overlap was found in B. Restricted by -f and -r.


Does it work without a space between the ``` and the {r}? It may be that the Rmd processor is getting confused about the spaces, but ignoring the bare r, and assuming the code is R by default.
The examples here don't have spaces:
https://rmarkdown.rstudio.com/lesson-3.html


The mroast function has an argument to specify which contrast do you want to test, quoting from the help page:

  contrast
  contrast for which the test is required. Can be an integer specifying a column of design, or the name of a column of design, or a numeric contrast vector of length equal to the number of columns of design.

So if you have the contrast of 1-min timepoint positive vs 1-min timepoint negative you can use this.
Read the help pages in detail, specially the help pages in limma are quite extensive and detailed to answer most questions if not all about how to use the software.


When running ChIPQC bioconductor package after peak calling, should I use sorted BAM files or unsorted BAM files?

You need to provide sorted bam files, ideally with an index. If you run ChIPQCsample, it calls ChIPQC:::sampleQC which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.
We can use an example, first in linux:
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/wgEncodeHaibTfbsA549RxlchV0416101Etoh02AlnRep1.bam -O Etoh02AlnRep1.bam

samtools index ./Etoh02AlnRep1.bam
# s
samtools view -uh Etoh02AlnRep1.bam chr1 > sorted.bam
samtools index sorted.bam
# quick way to unsort
samtools sort -n -o unsorted.bam sorted.bam
# it fails to index
samtools index unsorted.bam
[E::hts_idx_push] Unsorted positions on sequence #1: 246814554 followed by 42922635
samtools index: failed to create index for "unsorted.bam"

Then in R:
library(ChIPQC)
ChIPQCsample("unsorted.bam")
Bam file has 25 contigs
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: unsorted.bam

ChIPQCsample("sorted.bam",annotation=NULL)
Bam file has 25 contigs
Calculating coverage histogram for chr1

Calculating SSD for chr1

Calculating unique positions per strand for chr1

Calculating shift for chr1

 300 / 300
Counting reads in features for chr1

....


Standalone Blast+; Automating searches, Blastn formatting, interpretation

You just have to change the output format using the -outfmt parameter. To have an output in a tabular format just use -outfmt 6. The resulting table have 12 columns containing the following informations queryacc subjectacc pident length mismatch gapopen querystart queryend subjectstart subjectend evalue bitscore. You may restrict the number of hits per query using -max_target_seqs 5 for example.


Summing two columns where they share interval overlap

If you want to sum two columns based on overlaps, converting them to RLEs is a good idea.

# get the cpg scores only where there is overlap with the genes you are interested in
cr1 = cpg1.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")
cr2 = cpg2.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")

# convert these to run length encodings, using the score from the CPG column
rle1 = cr1.to_rle("CPG", rpm=False)
rle2 = cr2.to_rle("CPG", rpm=False)

rle_sum = rle1 + rle2

# convert the rles back to ranges
# (I will give the to_ranges function an option dtype eventually)
range_sum = rle_sum.to_ranges().apply(lambda df: df.astype({"Start": np.int32, "End": np.int32}))

# add the gene names to the scores
gr = range_sum.join(regulatory_elements) # .drop(like="_b")

Now the result is:

print(gr)
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 15623     | 15723     | 77.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 16034     | 16134     | 84.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 36588     | 36688     | 2.0         | +          | 36479     | 46479     | +            | 2871      |
# | chr1         | 37271     | 37371     | 30.0        | +          | 36479     | 46479     | +            | 2871      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59312851  | 59312951  | 28.0        | -          | 59308107  | 59318107  | +            | 98599     |
# | chrY         | 59332843  | 59332943  | 82.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59337592  | 59337692  | 41.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59341719  | 59341819  | 15.0        | -          | 59332664  | 59342664  | -            | 99215     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 668,317 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

print(gr[gr.Score > 100])
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 96874     | 96974     | 116.0       | +          | 87668     | 97668     | -            | 7391      |
# | chr1         | 96874     | 96974     | 116.0       | +          | 89794     | 99794     | -            | 6357      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 744041    | 754041    | +            | 1150      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 749277    | 759277    | -            | 6395      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59052942  | 59062942  | -            | 99957     |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59063570  | 59063576  | 105.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59269010  | 59269085  | 152.0       | -          | 59264704  | 59274704  | -            | 99499     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 106,643 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.


Does this FASTQ data contain single or paired end calls?

Entries in a fastq file occupy 4 lines each and R1 and R2 are typically in different files. Since that SRA project is 2x51 you seem to have run fastq-dump without the --split-3 option, so both R1 and R2 are merged together. Make your life easier and never use SRA, but instead just download the individual files from ENA.


This is a common question, and the answer is while you can calculate this, it's not statistically robust.
You're likely to arrive at false conclusions.
Why is this?
Because while TPM is something you can calculate from a single sample, your question involves multiple replicates and those replicates vary between each other.
Counts in units of TPM are missing information about the variation between samples.
This extra information affects fold change estimates, and is why differential analysis of sequencing data is complicated.
See this answer to a related question, and this often-cited blog post for more information.


The gi is an abbreviation for "Genbank identifier".
This is a pretty standard convention used by data stored in NCBI databases. There used to be a pretty comprehensive description of the conventions used at NCBI (I wouldn't say it was a standard or specification, just convention) here, but this page is no longer available it seems.
The AY848686 label is an accession value, and the appended number is a version number.

See this page for more info.

The short answer is: these are all internal conventions. None of it needs to be there, except perhaps to ensure NCBI's internal tools work as expected.

UPDATE: I found a non-official copy of the old docs here. Look for the section starting with "FASTA Defline Format". I've copied it here to prevent further problems with bitrot in the future.

FASTA Defline Format
--------------------
The syntax of FASTA Deflines used by the NCBI BLAST server depends on
the database from which each sequence was obtained.  The table below lists
the identifiers for the databases from which the sequences were derived.

  Database Name                         Identifier Syntax

  GenBank                               gb|accession|locus
  EMBL Data Library                     emb|accession|locus
  DDBJ, DNA Database of Japan           dbj|accession|locus
  NBRF PIR                              pir||entry
  Protein Research Foundation           prf||name
  SWISS-PROT                            sp|accession|entry name
  Brookhaven Protein Data Bank          pdb|entry|chain
  Patents                               pat|country|number
  GenInfo Backbone Id                   bbs|number
  General database identifier           gnl|database|identifier
  NCBI Reference Sequence               ref|accession|locus
  Local Sequence identifier             lcl|identifier

For example, an identifier might be "gb|M73307|AGMA13GT", where the "gb" tag
indicates that the identifier refers to a GenBank sequence, "M73307" is its
GenBank ACCESSION, and "AGMA13GT" is the GenBank LOCUS.

"gi" identifiers are being assigned by NCBI for all sequences contained
within NCBI's sequence databases.  The 'gi' identifier provides a uniform
and stable naming convention whereby a specific sequence is assigned
its unique gi identifier.  If a nucleotide or protein sequence changes,
however, a new gi identifier is assigned, even if the accession number
of the record remains unchanged. Thus gi identifiers provide a mechanism
for identifying the exact sequence that was used or retrieved in a
given search.

For searches of the nr protein database where the sequences are derived
from conceptual translations of sequences from the nucleotide databases
the following syntax is used:

                     gi|gi_identifier

An example would be:

        gi|451623           (U04987) env [Simian immunodeficiency..."

where '451623' is the gi identifier and the 'U04987' is the accession
number of the nucleotide sequence from which it was derived.

Users are encouraged to use the '-I' option for Blast output which will
produce a header line with the gi identifier concatenated with the database
identifier of the database from which it was derived, for example, from a
nucleotide database:

        gi|176485|gb|M73307|AGMA13GT

And similarly for protein databases:

        gi|129295|sp|P01013|OVAX_CHICK

The gnl ('general') identifier allows databases not on the above list to
be identified with the same syntax.  An example here is the PID identifier:

        gnl|PID|e1632

PID stands for Protein-ID; the "e" (in e1632) indicates that this ID
was issued by EMBL.  As mentioned above, use of the "-I" option
produces the NCBI gi (in addition to the PID) which users can also
retrieve on.

The bar ("|") separates different fields as listed in the above table.
In some cases a field is left empty, even though the original
specification called for including this field.  To make
these identifiers backwards-compatiable for older parsers
the empty field is denoted by an additional bar ("||").


Duplicate gene symbol handling in GEO gene expression data

Since this is a microarray (U133A), what you're seeing is that each probe is associated with one or more transcripts. The general strategy when dealing with microarrays is to use RMA (or fRMA) to summarize to probesets and then collapseRows to summarize to gene-level information. Since this is all typically done in R, you can find some further discussion on the bioconductor site. I would suggest going through this process and then only using the result of that in building your network (whether you're happy with sticking to transcript-level quantifications or gene-level quantifications for building your network is up to you).


How do I download the mitochondrial haplogroup datasets for human genetics online?

Hmmm, I am not experienced enough in biology or bioinformatics to truly understand the question, but I am guessing you want to download the FASTA file or genebank file containing the sequence of the accession number you listed above.
There are 2 ways to do this -

Using your browser, if you want to download a small sample of files. On the right hand side of your screen you can see the "send to" button on this link, click on it and choose the what's appropriate for you and then click on create file.

https://www.ncbi.nlm.nih.gov/nuccore/MT742594.1

If you familiar with Biopython and python programming you can download multiple files at once in whatever format you want, there is an easy script to do this and has been taken from Biopython's documentation on biopython.org
import os
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are
filename = "MT742594.1.gb"
if not os.path.isfile(filename):
 # Downloading...

 net_handle = Entrez.efetch(

     db="nucleotide", id="MG762674", rettype="gb", retmode="text"

 )
 out_handle = open(filename, "w")

 out_handle.write(net_handle.read())

 out_handle.close()

 net_handle.close()

 print("Saved")

print("Parsing...")
record = SeqIO.read(filename, "gb")
print(record)

You can make a list of the accession numbers of the sequences you want to download and pass the list in the id, if you are curious as to how the code works visit biopython.org .


Calculating p-value for introns which are retained (or not) between 2 conditions

Here is how we do this:

First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.

transcript 1 |>>>>>>>>>|----------|>>>>>>>>|---|>>>>|
transcript 2 |>>>>>>|-------------|>>>>>>>>|---|>>>>|
transcript 3 |>>>>>>|-------------|>>>>>>>>>>>>>>>>>|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|

this can be achieved using cgat gtf2gtf --method=genes-to-unique-chunks from the cgat-apps toolkit if you don't want to code it yourself.

I call
chunks 1, 4 and 6 constitutive exon
chunk  3 constitutive intron
chunk  2 alternate intron
chunk  5 annotated retained intron.

Now use featureCounts to count the number of reads in each chunk.

Test the chunk counts for usage using DEXSeq.

Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.


How can I download genePred file from UCSC?

It turned out that using table browser on UCSC we could download an output with all fields. That output actually contains a GenePred file but has one extra column on the left end.

The genePred file format definition is here.


Changing this code in a way to work for my files

Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible.

Here's a Perl script that should do exactly what your bash script does but will be several times faster and can deal with multiple input files:

#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile) or
      die("Couldn't open '$outFile' for writing: $!");
  open(my $inFileHandle, '<', $inFile) or
      die("Couldn't open '$inFile' for writing: $!");

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    my @infoFields = split(/;/, $info);
    my $readCount = $infoFields[12];
    $readCount =~ s/ReadCount=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCount = $infoFields[25];
    $variantAlleleCount =~ s/VariantAlleleCount=//;
    #/#Useless comment to fix SE syntax highlighting
    my $referenceAlleleCount= $readCount - $variantAlleleCount;

    my $readCountControl = $infoFields[13];
    $readCountControl =~ s/ReadCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCountControl = $infoFields[26];
    $variantAlleleCountControl =~ s/VariantAlleleCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $referenceAlleleCountControl = $readCountControl - $variantAlleleCountControl;
    my $VAF = $infoFields[27];
    $VAF =~ s/VariantAlleleFrequency=//;
    #/#Useless comment to fix SE syntax highlighting

    my $outString = join("\t", $outFileName, $chrom, $pos, $ref, $alt, $readCount,
                         $variantAlleleCount, $referenceAlleleCount, $readCountControl,
                         $variantAlleleCountControl, $referenceAlleleCountControl, $VAF);
    print $outFileHandle "$outString\n";
  }
  close($outFileHandle);
  close($inFileHandle);
}

Save that file as parser.pl and then run it on all your vcf files with:

perl parser.pl /path/to/vcfs/*vcf

That should create a .parsed.txt file for each input file in the target directory.

If you just have to use the shell for this, you can try with this improved version of your script:

#!/usr/bin/env bash

for file in "$@"; do
  outFile="${file}.parsed.txt"
  outFileName=$(basename "$outFile")
  (
    printf '%s\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t' "$outFileName"
    printf 'TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount\t'
    printf 'NormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n'
  ) > "$outFile"

  echo "Reading $file"
  grep -v '^#' "$file" |
    while read -r line; do
      read chrom Pos Ref Alt rest < <(awk -F"\t" '{print $1,$2,$4,$5}' <<<"$line")
      ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
      VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
      ReferenceAlleleCount=$(($ReadCount - $VariantAlleleCount))

      ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
      VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

      ReferenceAlleleCountControl=$(($ReadCountControl - $VariantAlleleCountControl))

      VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
      printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$outFileName" "$chrom" "$Pos" "$Ref" "$Alt" \
             "$ReadCount" "$VariantAlleleCount" "$ReferenceAlleleCount" "$ReadCountControl" "$VariantAlleleCountControl" \
             "$ReferenceAlleleCountControl" "$VAF" >> "$outFile"
    done
done

Just to illustrate what I mean about the shell being slow, I ran both scripts above on the example files you provided:

$ time foo.pl *vcf 2>/dev/null

real    0m0.051s
user    0m0.047s
sys     0m0.003s

$ time bar.sh *vcf 2>/dev/null

real    2m28.842s
user    3m33.648s
sys     1m12.456s

As you can see above, the perl script took less than a second, while the bash script took almost two and a half minutes!

Finally, if you find the perl script too long, here's the same thing in a more condensed style:

#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile);
  open(my $inFileHandle, '<', $inFile);

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    $info=~s/[^;]+=//g;
    my @infoFields = split(/;/, $info);
    print $outFileHandle (join "\t", $outFileName, $chrom, $pos, $ref, $alt,
                               $infoFields[12],$infoFields[25],
                               $infoFields[12]-$infoFields[25],
                               $infoFields[13],$infoFields[26],
                               $infoFields[13]-$infoFields[26],
                               $infoFields[27]) . "\n";
  }
}


I think there is no general answer to your question.

Which type of ID you will use depends on how practical it will be for your downstream processing. This depends on what other source of information you have, how much conversion work you will need to do in order to match your data to these information if you use one or the other.

I think the choice may emerge by trial and error, as you realize that a given choice is actually not convenient in your particular use case. Try to set up easy to use conversion procedures that you can easily insert in some steps of your workflow, so that you can easily switch from one choice to the other.

Concerning the tools to do the conversion, see the following question and its answers: Converting gene names from one public database format to another


Seurat FindMarkers() output interpretation

The p-values are not very very significant, so the adj. p-value. You need to plot the gene counts and see why it is the case. It could be because they are captured/expressed only in very very few cells. VlnPlot or FeaturePlot functions should help. You haven't shown the TSNE/UMAP plots of the two clusters, so its hard to comment more.


Where to download baseline/average gene expression level of all human coding genes?

There is not going to be a public database with this data. Apart from anything else, generating data for early human development is difficult and ethically tricky. Also, most people who care about development would probably think that a whole body average would be meaningless.

The best you could do would be to select tissues you were interested in from something like Expression Atlas and calculate the averages yourself.


how to set database other than nr for remote blast+ search

I believe you're looking for env_nr? It's listed as such, under Metagenomic proteins in the blastp webpage. It appears that the word within brackets should be supplied alongside the -db parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession.

Edit: I followed up on this with a little more digging around, and it appears the database you're looking for is called Microbial_proteins. I ended up quote-searching Microbial proteins from nr on Google and exactly 2 hits turn up. The first hit, I think, is what the OP posted a screenshot of in a comment to this answer. The second is a list of databases along with keywords from a Spain-based bioinformatics firm: http://data.biobam.com/ncbi_blast_dbs_protein.pdf. I retried the dummy example with the Microbial_proteins database, and it seems to work. Again, I am not sure why I can't easily find this on an NCBI resource. Also, I am curious as to how different this database is from env_nr.


Protein PDB/CIF text file reader that can fold a same section into one line

Programs that work with PDB and mmCIF files may present the structure as a tree, for example:

(this is a Find dialog from ccp4MG)

but I'd be surprised if a text editor had such a feature. Normally, you don't edit PDB files in a text editor.


Difference between genome assembly and genome sequence alignment to a reference to find structural variants

the first being what the 1KGP and SDGP did and sequence the whole
genome, align the reads to the reference genome and end up with a BAM
file.

If you have a well defined reference genome (e.g. human, mouse etc.) and you are interested in population level genetic variation, then this is the main approach. If you sequence a new human genome in the classical way (i.e. short reads ~30X coverage etc.), de novo assembly is pointless generally and you can more rapidly map reads to the reference. The read mapping approach has the large advantage that when you get the to variant calling stage, you can use information about both the depth of sequencing and the base quality scores.

The second is to assemble personal genomes and then compare or align
the assemblies to each other and the reference genome or using the
Lastz/LiftOver/ChainNets

This is the traditional comparative genomics approach used for non-model organisms and for conducting evolutionary rate comparisons at the cross-species level. You discard information on sequencing depth and quality scores, but you can compare many species' genomes at once genome-wide. You also worry less about sequencing errors etc. because you are not looking for relatively rare population level variation, but instead more common (depending on your species) variation between species.


What are doublets in single cell RNA-seq data?

Is doublet a set of cells sequenced as a single cell?

Yes. Depending on the method of single cell sequencing it may be more or less likely for groups of cells to be captured and barcoded with the same "unique" barcode. This is more likely in split-pool RNA sequencing (e.g. SPLiT-seq), and less likely in cell-capture RNA sequencing (e.g. Fluidigm C1). Doublets can also be created through physical / experimental processes (e.g. from tissues that were not completely dissociated).
Bear in mind that detecting doublets is not as simple as counting for doubling of transcript expression, because the expression profiles are different in different cells. That's why single-cell sequencing is useful in the first place.


I'm not sure what you mean by "trim". The code below will discard any reads that have a 10bp window with a mean phred score below your cutoff. I'll ammend the post to also include a solution that "trims" basepairs on both ends until the window achieves an average score of above the cutoff. It would be ideal if you could provide a minimal reproducible example that shows what your code does and what you'd like it to do.
As far as the code goes, I'd first suggest using a library like Biopython to parse the Fastq file. It will make your life a lot easier. For your code, it looks like seq actually contains all the sequences, so len(seq) is the number of reads in the file. However, when you iterate over i in range(len(seq)), it looks like you're using i to both index which read to iterate over as well as get the 10 base pair window. You need an additional index variable to obtain the window:
for window_start in range(len(seq[i]) - window_len) + 1)
    window = seq[i][window_start:window_start + window_len]

Additionally, some pythonic tips for the future:

If a range starts from 0, don't specify it as an arugment i.e. use range(x) as opposed to range(0, x)
Try to avoid anything along the lines of range(len(x)). In this case, you could just iterate over the sequences directly: for s in seq: or, if you wanted the index, you could enumerate the sequences via for seq_idx, s in enumerate(seq)
Use good variable names. seq implies a single sequence, yet is a list of sequences. Same with qual

Here is some code that will retain only sequences that have all windows with a passing mean phred score
from Bio import SeqIO

read_file = "/file/here.fq"
w = 10
cutoff = 15
valid_seqs = []
with open(fq_file) as fq_fd:
    seqs = SeqIO.parse(fq_fd, "fastq")
    for seq in seqs:
        scores = seq.letter_annotations["phred_quality"]
        mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
        if any(mean_score < cutoff for mean_score in mean_scores):
            print(f"{seq.id} does not pass the filter!")
            continue
        valid_seqs.append(seq)


How to extract reads with INDELs > a given size?

There's no simple way to do that with awk, since you would need to parse the CIGAR string and iterate over the operations. In python using pysam you would iterate over the reads and filter as follows:
keep = False
for OP, LEN in b.cigartuples:
    if OP == 2 and LEN > 20:
        keep = True
        break

if keep:
    ...do stuff...


You could use Biopython. If its a fasta file, it will be complicated to write fasta output (especially multiline fasta) with grep or awk. Simple solution is to use biopython, so that you can even match for any complex patterns in the fasta header.
from Bio import SeqIO

rna_records  = []
dna_records = []

for seq in SeqIO.parse("in.fa","fasta"):

    if "RNA" in seq.id:
        rna_records.append(seq)

    elif "DNA" in seq.id:
        dna_records.append(seq)

SeqIO.write(rna_records, "RNA_out.fa","fasta")
SeqIO.write(dna_records, "DNA_out.fa","fasta")

This code needs to be optimised if you are dealing with large sequences like human genome fasta as it holds everything in memory.


Is this a low complexity region in our human genome?

Yes, that's a short low complexity region wedged between a SINE and an snRNA. More importantly, your alignments have MAPQ of 0 (that's why they're filled with white in IGV), which will happen if they map equally well to multiple locations. Without looking at the sequence one can use that alone to determine that these are not trustworthy mappings.


Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?

The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power.

As this question suggests, if only a specific field from alignments is needed (or fields), then it will probably be better to extract just those fields and do any necessary conversion only on those fields. While this can be done by writing a BAM parser from scratch, many developers have already written software libraries to process BAM files in this way.

Devon Ryan has suggested htslib, which is a C library written by a group that includes the people who wrote the SAM/BAM/CRAM file format specifications. There's also pysam, which is a python wrapper around htslib.

The particular tool that is used will depend on your familiarity with programming and the specific thing that you want to do. If you just want to "output a BAM file", then cat is the quickest:

cat file1.bam

... but you probably don't want to do that, because it seems like you want to process a text-based representation with a script. Because you haven't specified which fields you are interested in, it's not possible to suggest the best thing to use. In the end, I expect that an "efficient" solution to your problem would involve htslib in some form. This is not really hacking samtools, it's using the backend of samtools to process BAM data.

However, it's useful to also think about the time cost of coding. What application do you have which means that the text-based processing method is not sufficient? It takes time to write code, and a lot more time to debug that code to make sure it's doing the right thing in all situations. If this is for a one-off thing, then samtools view output fed into your script may be the quickest solution.


How to interpret contig-alignment.psa produced by velvet

It looks like you are trying to do a transcript isoform expression analysis. Velvet is a genome assembler, and really not designed for this task. The tools you are looking for are kallisto and salmon.


How can I do an overlapping sequence count in Biopython?

For Biopython 1.70, there is a new Seq.count_overlap() method, which includes optional start and end arguments:

>>> from Bio.Seq import Seq
>>> Seq('AAAA').count_overlap('AA')
3
>>> Seq('AAAA').count_overlap('AA', 1, 4)
2

This method is also implemented for the MutableSeq and UnknownSeq classes:

>>> from Bio.Seq import MutableSeq, UnknownSeq
>>> MutableSeq('AAAA').count_overlap('AA')
3
>>> UnknownSeq(4, character='A').count_overlap('AA')
3

Disclaimer: I co-contributed the .count_overlap() methods with Peter Cock, see 97709cc


Integration of different microarray dataset to run GSEA

It makes more sense to evaluate by separate the pathway or gene set you want and see if in the three datasets result in a coherent message than to merge these datasets, as you will mix batch effects and treatment effects.

When in several datasets one looks for the same statistic it is usually considered meta analysis, not integration.


I think this is more to do with tree building than alignment construction.

I believe the underlying issue is the 1) distribution of indels, this results in a blockwise deletion in the resulting phylogeny and 2) nucleotide saturation.

Indels Thus, for indels

Archaeal
XXXX----YYYYY
bacterial
----ZZZZYYYYY

The resulting phylogeny will be based on YYYYY and XXXX is ignored.

The alternative explanation is that the genetic distances involved have resulted in saturation of the nucleotide sequences. SSU has to be a nucleotide phylogeny and the level of reversion mutations over a vast amount of evolutionary time can become indistinguishable from random.

With regards the first point, you can replace the indels "-" with "N". This is a cheat to trick a tree building program into not discarding the indel. Model based tree building programs have no ability to create a genetic distance between a base pair and an indel, so this data is essentially discarded from the tree. Worse it will remove the site from phylogenetic inference - which can be a massive loss of data. Two notes one this,

Some max. likelihood, Bayesian tree programs can avoid blockwise deletion, I don't know programmatically how they do this;
In PAML the distribution of indels was modelled, it doesn't appear to have been widely adopted however, but might be of value.

Alternative The alternative method is really for protein genes and called "RY" coding, it might work on SSU because of the compensatory mutations involved in maintaining a secondary structure of the SSU.

Word of warning, the saturation may be so great, very little nucleotide point mutation is preserved. Essentially we are talking about preservation of phylogenetic signal across the origins of life ... and that is a very, very long time ago.


FASTQC overrepresented sequences after trimming

As @AaronBerlin mentioned, you didn't remove reads that were completely trimmed. Next time use the --minimum-length option and set it to something reasonable, like 20. Alternatively, use "Trim Galore!", which is a wrapper around cutadapt that has more reasonable defaults.


Seurat Violin Plot: Why do dots align in one row?

That is because SCTransform (mostly) outputs discrete values, you can check the @slot of your Seurat object.


residual Squared Coefficient of Variation (rCV²) vs Distance to Median (DM)

There has been a lot of work done on this problem already, particularly these two papers:

https://www.biorxiv.org/content/10.1101/576827v2

https://www.biorxiv.org/content/10.1101/574574v1


What value of reference RMSD after docking simulation from any starting coordinate would be considered outrageous to have?

Let's get some definitions out of the way:

RMSD, root mean square deviation, is a metric of distance of molecule A and molecule B. Think of it as average euclidean distance. It us casually called just 'deviation', but has nothing to do with standard deviation (the root of variance, 2nd moment). The units are Ångström.

Reference RMSD is a term that appears in some manuals and is a the RMSD between the calculated pose and the starting pose. When you set up your docking you as a user define the starting position and how much it can deviate from that position. If you have a crystal structure with a ligand and you are docking the same or similar compound, then yes, it makes total sense to worry about this metric. In which case, RMSD greater than 2 are worrisome in a 2 Å resolution structure. You especially see RMSD as a metric of accuracy when PDB structures are redocked, wherein the same ligand is docked and reference is the file from the PDB (or its PDB Redo derivative), so RMSD greater than 2 Å is bad as it would not match the density, while under 1 is really good —note that in tests about 50% and 75% of re-docked ligands are docked less than 2 with AutoDock (free) and ICM-Dock (pay) respectively. But if you are doing it for an experiment and not for a test, you would probably have set the docking to be highly constrained in the first place (the grid size in AutoDock or position randomisation factor in most others, eg. Rosetta Ligand dock). If you actually care about reference RMSD the thing to manually look at is the top 3-4 poses in the electron density of your native structure. If you don't really know how your ligand docks please ignore this value as it is just a product of your starting pose. It is not a measure of how "happy" the complex is.

ΔGibbs is a negative value in kcal/mol that is basically a metric of how "happy" the atoms are. It is actually an energy potential, but confer Wikipedia pages for the proper definition. The more negative it is the happier the atoms are. Force-field calculations solve this value and the difference between undocked and docked (ΔΔG) is a good metric of how well your ligand docks. Please use this.

Affinity ΔΔG is not a k_off or a k_on. That is a result of the ligand binding to nearby positions and how well it can get into the binding site (e.g. in P450 enzymes the ligand orientation is dictated by how it squeezes into the entrance way). And how insoluble the ligand is, which is a major problem*. For a better metric, there is dynamic undocking where a constrained MD simulation is done with the ligand getting pulled out of the pocket. But this is not trivial.

∗ As any comp-chemist can tell you it is mighty hard to calculate solubility. Some force-field calculations do not even take account of π-π or S-π interactions. Drug-like compounds are rigid aromatic and rigid and highly hydrophobic. Many programs use a hack called implicit waters, where the water is not simulated as individual molecules but as a flowing thing. This means that ΔΔG benefits of constrained waters are ignored and alternative protontation states of water too. This is fine for protein (mostly), but has major errors for drugs. But in terms of thing to look out for in a docking algorithm, explicit waters is not the first as conformer sampling and side chain perturbation are more important.


Is there a aminoacid mutation that forms similar protein?

A mutation is most often likely to be either neutral (no effect) or destabiling (protein misfolds). Changes in activity require changes at the active site or the entranceway to it —which are often only a handful of residues, so are very uncommon.

So the answer to your question is most mutations keep the activity, although a large fraction break it. While to get a new activity the process is more fluid and is not boolean.

NB. Not all proteins are enzymes, but for simplicity I will concentrate on enzymes.

Neutrality and core vs. surface

The outer residues of a protein, the surface residues, are very tolerant to mutations. Unless they are part of an binding surface (interface) or the change results in a charge change that is. The core residues on the other hand are more sensitive to mutations.

Destabilisation

One can do a scan in silico or in vitro or in vivo to determine the mutational landscape of a protein, which can take various forms.

This being Bioinformatics stack overflow I'll give an in silico example. The following for example is a heatmap I did some time back with Rosetta pmut_scan, which is a classic tool for this and shows in the y axis the residue it was mutated to and on the x axis the position. A larger change in Gibbs free energy of folding) relative to wild type (∆∆G) means it is destabiling, but anything less than 5 kcal/mol is neutral. As you can see many are neutral, but some half are destabiling —with a caveat for proline, which require a backbone change.

In the in vitro setting one speaks of melting temperature (DSC assays) or turnover number over Michaelis constant (enzyme specific assay). In the in vivo setting one speaks of a fitness cost of a mutation, i.e. the protein is less active so the organism grows less.

Drug resistance

Drugs inhibit enzymes by binding the active site. Therefore resistance is often costly. There are some nice papers that look into the fitness cost and mutational landscape of HIV protease or β-lactamase. In this situation a strong epistasis is seen. That is, a first mutation that confers resistance to the drug greatly destabilises the protein, while a second compensates for the first mutation.

Change in function

An enzyme can catalyse a reaction well on a given substrate, because it was evolved this way (technically the encoding gene evolved, but shh). This is the physiological activity. However, enzymes can also catalyse side activities on similar substrates with a lower turnover and higher Michaelis constant. This is especially true for substrates it was not selected against. This is called "enzyme promiscuity". The tolerance of the physiological activity to mutation is called "robustness", while the ability to increase promiscuous activities is "plasticity". Different enzymes have a different balance of the two.

So a mutation may change the strength of the physiological activity, but may result in a new promiscuous activity that can be used as a stepping stone for selection to increase this activity. Therefore, as you can see it is not a binary switch, except for a few very specific case.

Best amino acid for changes?

Proline is a very strange amino acid as it is a secondary amino acid. It lacks a hydrogen on the backbone nitrogen so it cannot form the required hydrogen bond for an α-helix. And can adopt a cis-conformation with minimal penalty (whereas it is about 5 kcal/mol for other chiral amino acids). So changes to this amino acid are often very destabilising.

Protein cores tend to be hydrophobic so changes to the hydrophobic residues are often deleterious, especially when going from a small residue to a larger one.

The remaining residues are frequently found on the surface so may be neutral. However, some are found in the active site and may be able to shift the parameters of the ligand space for that enzyme. A change in charge may result in an increased specificity for a different set of substrates for example. A size change from a larger amino acid to a smaller one may result in a broader range of promiscuous activities —mutations to alanine for example.

So ironically, polar and charged residues as they are frequently found on surfaces are the residues which are most often neutral, but in a small amount of cases they may change activity.


You have scale="col" in the code. What you are plotting is the z-score calculated based on the value distribution per column. Try changing it to "row" or "none".


How to reduce the occupied RAM when you are dealing with a very sparse matrix in a single-cell Experiment in R?

Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well.
Anyway, the other answers and comments are way off. scran has supported sparse matrices for years, ever since we switched over to the SingleCellExperiment class as our basic data structure. quickCluster does no coercion to dense format unless you tell it so explicitly, e.g., with use.ranks=TRUE (in which case you're asking for ranks, so there's little choice but to collapse to a dense matrix).
You don't provide an MWE or your session information, but this is how it rolls for me:
# Using the raw counts in the linked dataset. Despite being
# called a CSV, it's actually space delimited... typical.
library(scater)
mat <- readSparseCounts("GBM_raw_gene_counts.csv", sep=" ")

# Making an SCE just for fun. Not strictly necessary for
# this example, but you'll find it useful later.
sce <- SingleCellExperiment(list(counts=mat))

library(scran)
system.time(clust <- quickCluster(sce))
##   user  system elapsed
##  3.170   0.174   3.411

This is running on my laptop - 16 GB RAM but I'm definitely not using all of it. I only go full throttle when I'm working on some real data, e.g., the 300k HCA bone marrow dataset. Check out the book for more details.
Session info below, I don't know quite enough SO-fu to know collapse it.
R version 4.0.0 Patched (2020-04-27 r78316)
Platform: x86_64-apple-darwin17.7.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Users/luna/Software/R/R-4-0-branch/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/R-4-0-branch/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] scran_1.16.0                scater_1.16.1
 [3] ggplot2_3.3.2               SingleCellExperiment_1.10.1
 [5] SummarizedExperiment_1.18.1 DelayedArray_0.14.0
 [7] matrixStats_0.56.0          Biobase_2.48.0
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2
[11] IRanges_2.22.2              S4Vectors_0.26.1
[13] BiocGenerics_0.34.0

loaded via a namespace (and not attached):
 [1] beeswarm_0.2.3            statmod_1.4.34
 [3] tidyselect_1.1.0          locfit_1.5-9.4
 [5] purrr_0.3.4               BiocSingular_1.4.0
 [7] lattice_0.20-41           colorspace_1.4-1
 [9] vctrs_0.3.1               generics_0.0.2
[11] viridisLite_0.3.0         rlang_0.4.6
[13] pillar_1.4.4              glue_1.4.1
[15] withr_2.2.0               BiocParallel_1.22.0
[17] dqrng_0.2.1               GenomeInfoDbData_1.2.3
[19] lifecycle_0.2.0           zlibbioc_1.34.0
[21] munsell_0.5.0             gtable_0.3.0
[23] rsvd_1.0.3                vipor_0.4.5
[25] irlba_2.3.3               BiocNeighbors_1.6.0
[27] Rcpp_1.0.4.6              edgeR_3.30.3
[29] scales_1.1.1              limma_3.44.3
[31] XVector_0.28.0            gridExtra_2.3
[33] dplyr_1.0.0               grid_4.0.0
[35] tools_4.0.0               bitops_1.0-6
[37] magrittr_1.5              RCurl_1.98-1.2
[39] tibble_3.0.1              crayon_1.3.4
[41] pkgconfig_2.0.3           ellipsis_0.3.1
[43] Matrix_1.2-18             DelayedMatrixStats_1.10.0
[45] ggbeeswarm_0.6.0          viridis_0.5.1
[47] R6_2.4.1                  igraph_1.2.5
[49] compiler_4.0.0


Why does a missing label 11plex TMT shows up at almost 50% intensity compared to other labels?

The developers answered along these lines:
This might be a normalization artefact. Reporter is normalizing the ratios so that the median is 1. So, if the label is almost not present the noise gets increased. To check that you can turn of the normalization. This function is currently not documented. The developers are working on it though.

GitHub issue


Getting sequences of one transcript per gene out of annotation and genome

Here is the custom python script solution:

#!/usr/bin/env python3

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt->gene
# jg25690   jg25690.t1
# jg25690   jg25690.t2
# jg25691   jg25691.t1
transcripts2genes = dict()
gene2transcript_mapfilename = 'data/genome/transcripts2genes.map'
with open(gene2transcript_mapfilename, 'r') as g2t_file:
    for line in g2t_file:
        gene, trainscript = line.rstrip('\n').split('\t')
        transcripts2genes[trainscript] = gene

# construct a dictionary gene->list of transcripts
# this fasta file is output of default gffread -w
genes2transcript_seuqneces = defaultdict(list)
transcripts_fasta_filename = 'data/genome/transcripts.fasta'
transcripts_fasta_file = SeqIO.parse(transcripts_fasta_filename, 'fasta')
for seq_record in transcripts_fasta_file:
    gene = transcripts2genes[seq_record.name]
    genes2transcript_seuqneces[gene].append(seq_record)

genes_fasta_filename = 'data/genome/genes.fasta'
with open(genes_fasta_filename, 'w') as out_file:
    # iterate through genes
    for gene in genes2transcript_seuqneces.keys():
        # get all the associated transcripts
        transcripts = genes2transcript_seuqneces[gene]
        # find the longest, by default the one with index 0
        lengest = 0
        for index, tr in enumerate(transcripts):
            if len(transcripts[lengest].seq) > len(tr.seq):
                lengest = index
        # this would be nicer with Bio interface (with wrapping of sequence lines)
        # but I was lazy and blast does not care, only thing I do is converting sequence to uppercase (not even sure if necessary)
        fasta_string = '>' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)


I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.

It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was done or there was a bit of degradation. The primary thing to ensure is that the bias is similar across samples/groups. If that's the case then you don't have to worry that your downstream analyses are going be affected by some sort of per-group bias.


Anyway the chances are this is gene duplication. It will be a housekeeping type gene or else histone that sort of thing. You simply map the different genes onto the various metazoan classes. In your case nematodes and arthropods are basal metazoan, so their is correct given the clade they form part of all all representatives is that gene.

There is nothing weird you just need to read the deep phylogeneticists work like Tom Caveliar-Smith and Sandy Balfour, possibly Blair-Hedges

There are alot of eukaryotes prior to the invertebrates, just look at the Cambrian explosion


How to split FASTQ reads without re-running `fastq-dump`?

My comment above mentions the slow IO for the Python code, so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using zcat | awk with a switch-case on the line number to stream the data.

# first pass to generate _R1 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 1, 76); break; case 3: print "+"; break; case 0: print substr($1, 1, 76)} }' | gzip -nc > {output}_R1.fastq.gz

# second pass to generate _R2 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 77, 76); break; case 3: print "+"; break; case 0: print substr($1, 77, 76)} }' | gzip -nc > {output}_R2.fastq.gz

All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution.
If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.

But I will probably make use of fasterq-dump next time, and make sure I specify to split the read mates.


As you suspected, correlations will be heavily changed by such a transformation. The reason for this is also exactly as you described, namely that using ranks rather than values for the transformation eliminates all of the information needed to actually preserve these relationships. As an example, suppose we have a single observation with rank 1 across all observations that none the less shows a decrease in absolute value across them (i.e., there's a correlation with a negative slope over time). Such a transformation will result in a correlation of 0.

One might be able to find a different transformation that preserves this relationship, or might be able to simply pool all of the time points together before computing the ranks (this would then at least partially preserve some correlations, though depending on what's needed to properly normalize the data it may cause other problems), but in general there's no guarantee that this is possible.


How can the Autocovariances, autocorrelations, and autocorrelation coefficients be calculated from a Protein Amino Acid Sequence?

Protein autocorrelation is similar to usual autocorrelation of time series (that is, function $r(k)$ measuring correlation between values lagged by $k$), except that instead of time dimension we slide along the backbone of amino acid chain. The "dependent variable" can be any quantifiable property of the amino acids - then the autocorrelation measures similarity in that particular property.

See e.g. Fernandez et al., 2008, for more extensive background:
https://www.tandfonline.com/doi/full/10.1080/08927020802301920?scroll=top&needAccess=true