Generate ligands candidates based on protein shape

TL;DR: docking is much slower than any ML approach, but the ML approach can be constrained by pharmacophores dictated by the active site.

Side note: Scale

The scale for ligand space exploration is generally is generally several orders of magnitude higher than "hundreds": Zinc DB lists 750 million enumerated compounds, GDB13 has enumerate all possible compounds up to 13 atoms (970 million), Enamine Real is 1.2 billion etc. In an actual experiment, one would use a more restrictive set. As an example of a real project, in Covid19 moonshot project there are 2,800 user suggested follow-ups.

Side note: More than just binding

There are a few requirements for a drug to be effective, and binding is only one of them:

Must follow Lipinksi's rule of five, i.e. it must be membrane permeable and small.
Must bind.
Must be synthesisable. Whereas the experience of someone learnèd is by far the best away to tell, there are several machine leaning approaches that can filter a huge dataset. IBM actually have a fun online tool.
Must not have side effects

So pre-filtering the list of compounds, regardless of docked poses, will happened, be it manually, by cut-offs or by advanced machine learning.

CPU time: docking vs. machine learning

These various steps are generally followed in that order.
Calculating properties (MW, logP etc.) for a given SMILES string to filter based on Lipinski's rule takes a fraction of a second, so that is generally always a first step.

Whether the compound can synthesised or bought is a general left last and done manually. Pharmacology is complicated so is ignored.

Docking a ligand takes time and different approaches have different time scales. A decent algorithm that uses implicit solvent model takes about a minute per core. A short MD simulation takes an 1 hour per core or more. There are docking programs for screening that can be faster but these dock only one conformer of the small molecule against a rigid protein and the results are not of much value. Also setting up docking properly is not easy (cf. steps required for Autodock 4, each with their own caveats). So to recap what needs to be considered that results in the slowness:

Hundreds of conformers of the ligand
Some flexibility of the protein, side-chains repacking only or even backbone changes.
thousands of poses explored.
Optionally. Implicit water is not as good as a TIP3 water model

Machine learning have a much smaller time requirement and is way more attractive for publications. Generally these work by using ligand only information, such as properties, splitting it into pharmacophores, colour etc. This requires a decent dataset, such as empirical results, say from a fragment screen, what compounds bound and which didn't. Often there is a lot of hype with machine learning and the result is something close to something already seen, as well put in this blog post by Pat Walters. This is a side effect of the starting data. The recent Stokes et al. 2020 Cell paper uses AI perfectly because of their staggering automated empirical data collection pipeline.

ML after docking

However, ML can actually be used after docking to improve the quality of the score.
The Boyles, Deane, Morris 2019 Bioinformatics paper investigated ligand-only machine learning to see why do ML methods that rely solely on ligand somehow often work. It is a good paper from a sceptical computational biochemistry viewpoint.
The scorefunction used in docking by the sampler can be classified as forcefield-based, empirical or hybrid. In some cases empirical cases, many terms are ligand based properties (/features) that would normally go into a ML search (say one trained on failed hits and successful hits).
For example Autodock Vina (not 4) has regression derived factors in the final weights of the score

Structural features in machine learning: HotSpots

Information of the protein active site can be used in ranking compounds in a non-docking machine learning approach. Namely, finding out what kind of chemical groups bind preferentially to different parts of the protein's active site, i.e. making a "fragment hotspot map". Specifically, a hotspot map allows the determination of a pharmacophore, which is an abstract generalisation of a group of molecules with both both electronic and steric features —a coarse-grain model, if you will. This approach is also used to detect where the active site actually is in the extreme cases where it is not know (primarily a problem for pipelines, not humans). There are several implementations as described here, the first being FTMap (site, paper), which docks 16 different types of compounds. Waters can play a role and severely complicate matters, hence approaches can be taken to counter that. However, unlike docking protein rigidity is not a big nuisance: this is because a hotspot map is most often used for filtering/weighing by the "color" not the shape of the molecules. But shape based searching is possible.

Shape based scoring

In terms of shape-based searching, this is actually a question of having pre-computed conformers, which can be problematic (but can and is often done on the fly), and shape based scoring —for which there are different metrics, such as the classical RMSD to more complexly weighted scores such as SuCOS.


As you sequence a library you are sampling from the pool of all possible oligonucleotides found in it. This is a stochastic process, meaning that you are more likely to sequence a piece of DNA that is more prevalent.

As you increase the number of reads that you sequence for a library, you decrease the chances of seeing a new read that you have not yet seen before. If you plot the number of unique reads per total reads sequenced, it will look something like the cumulative density function of the exponential distribution, where you have diminishing returns on unique sequences as you continue to sequence.

If you sequence a library and quickly stop seeing unique reads, we say that this library has low complexity. Therefore, you have "oversampled" it if you sequence too deeply. I am guessing that your metagenomics library actually has a lot fewer species/unique sequences than they anticipated, and that you sequenced deeply anticipating more.

Edit: In response to your comments, the bigger problem in 16S analyses is not sequencing too deeply, it is over-amplifying your samples such that they no longer represent what is actually in your gut microbiome samples. In this case if you used read normalization you would not accurately represent what was in the library that you sequenced since the amplified molecules (16S) are so small that you only have len(16S molecule)-(read len) + 1 possible reads in each direction. You might consider removing PCR duplicates with GATK but this will still have the same effect of not representing what was in the sequencing pool. It is very possible that you have many PCR duplicates after 16S amplification. I generally don't trust amplicon sequencing for quantitive analyses, and don't really trust for low coverage "is-this-species-there-or-not" either, especially when dealing with the 'dark matter' of uncharacterized species. You might be well off just asking "is this species present in high coverage at this time point?"


coexnet package for network co-expression analysis, is looking for a filelist.txt that doesn't exist

Finally, I found it. filelist.txt file is accessible from the FTP server of the NCBI. For example for that specific microarray is stored here

So there might be a problem at software side and they don't get downloaded.


Which sequence alignment tools support codon alignment?

I don't know of any transcript-to-transcript aligners that are able to do this, but LAST can align transcript queries to protein reference sequences using a specified frameshift cost. Here's the specific documentation for that option:

  -F COST

  Align DNA queries to protein reference sequences, using the specified
  frameshift cost. A value of 15 seems to be reasonable. (As a special
  case, -F0 means DNA-versus-protein alignment without frameshifts,
  which is faster.) The output looks like this:

a score=108 s prot 2  40 + 649
FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD s dna  8 117 + 999
FFLQ-IKLWDP\STPH*IVSSP/PSDLISAHTLCPRMKSQDN

  The \ indicates a forward shift by one nucleotide, and the / indicates
  a reverse shift by one nucleotide. The * indicates a stop codon. The
  same alignment in tabular format looks like this:

108 prot 2 40 + 649 dna 8 117 + 999 4,1:0,6,0:1,10,0:-1,19

  The "-1" indicates the reverse frameshift.

I sent an email to the LAST mailing list about adding a frameshift penalty for transcript-to-transcript matching; I've been pleasantly surprised with the requested features that Martin Frith has added to LAST in the past. Unfortunately, in this case the problem is too difficult to sort out due to all the possible combinations that could happen, so it's unlikely to be implemented in LAST in the forseeable future (unless someone else writes that code).


How can I remove (non-trivial) duplicates from a VCF file?

It turns out that bcftools can do this (tested on bcftools-1.8), if you give it the reference genome to test against:

$ bcftools norm -d none -f hg19.fa test.vcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -d none -f hg19.fa test.vcf; Date=Wed Feb 27 16:08:44 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  .       C       T       .       PASS    .       GT      0/1
Lines   total/split/realigned/skipped:  2/0/1/0

For the more complex case of the multi-allelic variant in the second VCF example from the question, you can run it through bcftools twice. Once using norm to left-align and split multi-allelic variants, and then again to remove the duplicates:

$ bcftools norm -m -any -NO z -O v -o - ~/test2.vcf |
    bcftools norm -d none -f hg19.fa
Lines   total/split/realigned/skipped:  2/1/0/0
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -m -any -NO z -O v -o - test2.vcf; Date=Wed Feb 27 18:18:32 2019
##bcftools_normCommand=norm -d none -f hg19.fa -; Date=Wed Feb 27 18:18:32 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  529514  CAC     T       0       .       .       GT      0/1
chr12   529516  529514  C       T       0       .       .       GT      0/0
Lines   total/split/realigned/skipped:  3/0/2/0


How many reads do I need to cover the entire genome?

Let's see what your slides are claiming $\epsilon$ to be, using:

$$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$

We can rearrange to get

$$ \epsilon \approx \frac{G}{L} \times e^{-C} $$

Now, what does it mean for an entire genome to be covered? Having a clear understanding of this definition is crucial, and your misunderstanding seems to partially stem from this.

I claim it's natural to consider a genome as consisting of $G/L$ bins. In this setting, we have $N$ balls that each randomly goes into a bin. Then, the entire genome is covered when all bins contain at least $1$ ball/read.

From the Poisson distribution, we know that P(a bin has 0 balls) = $e^{-C}$.

Then, we define $\epsilon$ to be the probability that at least one of the $G/L$ bins has 0 balls. Given the way we've modeled the problem, this is exactly equal to

$$ \epsilon = 1 - (1 - e^{-C})^{G/L}. $$

From above, your slides are claiming $ \epsilon \approx \frac{G}{L} \times e^{-C}$. This might be strange, as it appears to be equal to the expected number of bins with 0 balls, instead of the probability that there is at least 1 bin with 0 balls. Furthermore, it's strange because it contains only one exponential, compared to our exact derivation which has a double exponential.

The mismatch in the number of exponential is a hint that another mathematical approximation is at play. We can use a first-order approximation of the binomial series to obtain

$$ (1-e^{-C})^{G/L} \approx 1 - (G/L)e^{-C} $$ which holds particularly when $|e^{-C}| \times \frac{G}{L} << 1$.

In your case, $\frac{G}{L} = 30,000$, so the approximation only really starts getting accurate when $C > 15$ or so.

Following through with our approximation, we get

$$ \epsilon \approx 1 - (1 - (G/L) e^{-C}) $$

$$ \epsilon \approx \frac{G}{L} \times e^{-C}. $$


How to calculate the number of cells in different cell clusters with Seurat v3?

The "identities" can be accessed with the Idents() function in Seurat v3.

You should change cell.num <- table(shox2crepos_mtmg@ident) to

cell.num <- table(Idents(shox2crepos_mtmg))


State of the art mutation simulation software

Question I suspect this question is about the phylogenetics of methylation and the approach the investigator is proposing would be the last approach to use.

Summary
The approaches to assess the phylogenetics of methylation in order of preference are:

dN/dS between CpG sites and non-CpG sites,
An explicit molecular clock between CpG sites and non-CpG sites
A likelihood ration test (LRT) based on a null distribution of mutations generated by a Monte Carlo algorithm

Background/Rationale
The approaches you are using have been around for a long while and generate the random distribution of mutations for a given phylogeny. Its phylogenetics in reverse and the technique is known as Monte Carlo, you start with a randomised probability and send it through a parameterised model to predict the amino acid/nucleotide. Thus it is ML(max likelihood), Bayesian, "phylogenetic-HMM" in reverse. It is used within likelihood ratio test calculations. The model is determined by a standard ML, Bayesian phylogeny algorithm, i.e. its circular mathematics because there is no independent means of calculating the mutation behaviour, so the precise context you are using this for needs to be carefully considered.

Packages There are lots of packages that do this type of Monte Carlo simulation, SOWHat being a good one (LRT) and you generate a large number of replicate data sets (100 or 1000). The one of the base algorithms is "seq-gen", although PAML may have implemented this.

Considerations To use this style of approach you need to carefully consider your question. The "mutation simulations" when analysed via a phylogenetic program will produce the same parameters that you initially set and trees of basically the same length. If you are using this to generate a null distribution for a phylogenetics test these approaches are useful, you then compare the observed likelihood against the null distribution. If you are using it to work out whether the mutation rate at CpG sites is higher than at other sites, its one approach amongst a number of alternatives.

Drawbacks Calculating a de novo null distribution of the mutation rate is computationally very expensive and therefore tends to the last calculation performed. It will address a single hypothesis, usually regarding topology.

Crux of the problem To really do what (I think) you want to do you need an independent measure of mutation and mutational behaviour and its not trivial to achieve. You would have to consider it like a machine learning calculation, with a formal training set, I don't know whether such an approach has been implemented.

Strengths/Summary In summary, Monte Carlo phylogenetics simulation has a limited application because of circularity of the calculation, BUT, BUT, BUT what I have omitted is when it is appropriate it is a very powerful test indeed.

Points of interest I do like the population dynamic simulations within Beast, this certainly "a priori" and its potential IMO needs further exploration. However, I don't think you looking at how molecular epidemiology might impact the mutational behaviour methylation. I have not looked at Indelible, but looks interesting.

I was going to talk about "processor war simulations" but its probably off-topic.


How to calculate fold change in gene expression from RSEM (RNAseq)

You should use a proper statistical framework for RNA-seq dfferential analysis (which includes FC calculation). Standard tools for this are (among others) edgeR or DESeq2. You could use tximport to import RSEM outputs into R and then use its output for e.g. DESeq2. The linked manual provides example code for this. Please read the manuals of edgeR or DESeq2 towards the respective downstream analysis. Both are extensive and compreshensive.


Add tags to the read pairs in a bam file

The best thing would be to figure out how to do this in UMItools. I don't know that tool at all, but it seems that it is supposed to be able to deal with / tag both reads in paired data. e.g. did you use --paired flag, etc.
I would suggest investigating that route first before you look at the next part of my answer.

Per Devon Ryan's comment suggestion I would suggest a workflow similar to this:

sort bam file by query name (samtools sort -n) to get the R1s immediately following the R2s
Write a short script/command that

reads through the bam file and for each line
looks up the tag in the R1, prints R1 line
reads the R2, substitutes the RX tag in, prints R2 line

pipe that through samtools view and possibly samtools sort back to a bam file of the same type that you started with.

My pysam is a little rusty but something like this might work
import pysam

bam = pysam.AlignmentFile("bam", "rb")
# print header
print(bam.header)

tag = None
for read in bam:
    r2 = False
    if tag:
        read.assign_umi_tag() # assign to read- not sure where/how
        r2 = True
    else:
        tag = read.lookup_umi_tag()  # i don't know where this tag is in bam format
    print(read)
    if r2:
        tag = None

I haven't tested that in any way, caveat emptor.
I am not personally familiar with UMItools and how it assigns tags so can't suggest anything more specific, sorry.


Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype?

Your approach is correct - it's not the only one permitted under their conditions, but it would be the standard way given only that information.

F-test compares a "full model" (predicting the phenotype using genotype and an intercept) versus a "reduced model" (predicting the phenotype using only the intercept, i.e. the mean value of phenotype). The p-value tells whether the full model predicts better, i.e. whether adding the genotype variable significantly improves the prediction. In this single-variable case, it should be approximately equal to the p-value for the beta coefficient of the genotype, if that's more familiar.

Anyway, this isn't specific to genetics, so you can check out other online sources - e.g. https://onlinecourses.science.psu.edu/stat501/node/295 for a longer, but accessible explanation of the F-test, and https://stats.stackexchange.com/ for many answers to similar questions.


How to separate a microbial strain of MTB resistant and susceptible drugwise?

Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl:

#!/usr/bin/perl

use strict;
use warnings;

## This will hold the information for each strain
my %strains;

## These are the functions that will test each strain
sub isMDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}){
    foreach my $drug (keys(%thisStrain)) {
      next if $drug eq "ISONIAZID";
      next if $drug eq "RIFAMPICIN";
      next if $drug eq "PYRAZINAMID";
      next if $drug eq "ETHAMBUTOL";
      next if $drug eq "STREPTOMYCIN";
      ## discard this strain if it is resistant to any other drug
      if ($thisStrain{$drug}) {
        return 0;
      }
    }

    return 1;
  }
  return 0;
}

sub isPreXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}) &&
  (
   $thisStrain{"CIPROFLOXACIN"}  ||
       $thisStrain{"MOXIFLOXACIN"} ||
   $thisStrain{"OFLOXACIN"}  ||
       $thisStrain{"LEVOFLOXACIN"}
      )
     ) {
    return 1;
  }
  return 0;
}

sub isXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  isPreXDR($_strain) &&
      (
       $thisStrain{"AMYCACIN"} ||
   $thisStrain{"CAPREOMYCIN"} ||
       $thisStrain{"KANAMYCIN"}
      )
     ){
    return 1;
  }
  return 0;
}

sub isSusceptible{
  my $_strain = shift;
  foreach my $drug (keys(%{$strains{$_strain}})) {
if ($strains{$_strain}{$drug}) {
      return 0
    }
  }
  return 1
}

## Iterate over each line in the input file.
## For each iteration, the special variable $_
# will hold the contents of the line.
while (<>) {
  ## Skip first line
  next if $. == 1;
  ## remove Windows-style line endings
  s/\r//;
  ## remove trailing newline
  chomp;
  ## split the line on commas
  my ($strain, $drug, $status) = split(/,/);
  if ($status eq "SUSCEPTIBLE"){
$status = 0
  }
  else {
$status = 1;
  }
  $strains{$strain}{$drug}=$status;
}

foreach my $strain (keys(%strains)) {
  print "$strain\t";
  if (isMDR($strain)) {
print "MDR\n";
  }
  elsif (isXDR($strain)) {
print "XDR\n"
  }
  elsif (isPreXDR($strain)) {
print "PRE-XDR\n"
  }
  elsif (isSusceptible($strain)) {
    print "SUSCEPTIBLE\n";
  }
  else {
    print "OTHER\n";
  }
}

If you save that script as foo.pl, you can run it on your file with:

perl foo.pl drug.csv

Selected output from running this on the drug.csv file you posted is:

$ foo.pl drug.csv | grep -E 'ERR038288|ERR038736|ERR038751|SRS553990|ERR133879'
ERR038751   OTHER
ERR038288   SUSCEPTIBLE
SRS553990   XDR
ERR038736   MDR
ERR133879   PRE-XDR

Various details about this question were added in chat: https://chat.stackexchange.com/transcript/message/47821646#47821646


Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value?

Let f and r be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base c in a proper 2-bit encoding, we update the two integers as follows:

f = (f<<2|c) & ((1ULL<<2*k) - 1)
r = r>>2 | (3ULL-c)<<2*(k-1)

With this updating rule, f keeps the forward strand k-mer ending at c and r is f's reverse complement. The hash of the k-mer can be min(f,r). If integer operations take constant time, computing the hashes of all k-mers of a sequence of length L takes O(L) time.

The above only works if k-mer fits a word, which means k<=32 on x86. For long k-mers, you can overload bit operators in C++. In case of 32<k<=64, a simpler solution is to split into lower and higher bits:

c1 = c&1, c2 = c>>1
f1 = (f1<<1|c1) & ((1ULL<<k) - 1)
f2 = (f2<<1|c2) & ((1ULL<<k) - 1)
r1 = r1>>1 | (1ULL-c1)<<(k-1)
r2 = r2>>1 | (1ULL-c2)<<(k-1)

This is twice as slow as the k<=32 version. It is possible to implement these lines with SSE2 intrinsics, but that is overkilling.

When you want to hash a long k-mer into fewer bits, there are a few options:

xor:      f^r
min:      min(f,r)
min+hash: hash(min(f,r))

where hash() is a generic randomization hash function, which could be Murmur3, Wang's integer hash function, etc.

The following is a microbenchmark. Here, we squeeze all k-mers into a 32-bit hash table. collision_rate = 1 - #distinct_hash/#distinct_kmer. This is not the best metric, but it somehow measures the randomness, which is better than nothing.

====================================================
 Algorithm    data  max-k  k   %collision  CPU time
----------------------------------------------------
 xor          chr11  32    31     8.5%       0.9s
 xor+Wang     chr11  32    31     9.6%       1.3s
 min+Wang     chr11  32    31     1.4%       1.3s
 min+Wang     chr11  64    31     1.4%       1.9s
 min+murmur3  chr11  inf   31     1.4%       5.7s
 min+Wang     chr11  64    51     1.5%       2.1s
 min+murmur3  chr11  inf   51     1.5%       6.8s
 min+Wang     chr11  64    63     1.5%       2.1s
 min+murmur3  chr11  inf   63     1.5%       7.5s
 radix sort   chr11  -     -      -        5.6-7.3s
----------------------------------------------------
 min+Wang     hg38   64    63    26.2%        37s
 min+murmur3  hg38   inf   63    26.2%       192s
 radix sort   hg38   -     -      -         ~138s
====================================================

Observations:

Don't use XOR. In fact, XOR hashes all palindromes to 0xffff.... This is already worrying enough.
Unless you want to work with very long k-mers, don't use generic string hash functions like FNV and Murmur. Murmur is even slower than radix sorting all hashes.
For k<=64, min+Wang is the best here. It is fast and simple to compute and has randomness nearly as good as murmur.


Pepsin digest (cleavage) does not work using RE?

I am really not sure if I've understood your rules correctly, but I think you want this:

If all of these are true:

position 2: not R, K or H
position 3: not P
position 4: not R
position 7: not P

And, if at least one of these is true:

position 1: P
position 1: P
position 4: L or W or Y or F
position 5: L or W or Y or F

Then cleave.

If so, then you probably don't want to do it using regular expressions. First because it will be hard, and second because it will be absolutely awful to maintain if you need to come back to it. Instead, try something like this:

#!/usr/bin/env python3
# X = input();
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
seq = list(X)
for i, res in enumerate(seq):
    if i > len(seq) - 7:
        break
    if (
            seq[i+1] in ('R','K','H') or
            seq[i+2] == 'P' or
            seq[i+3] == 'R' or
            seq[i+6] == 'P'
    ):
        continue
    if res == 'P':
        start = i
    elif seq[i+1] == 'P':
        start = i+1
    else:
        continue
    if (
            seq[i+4] in ('L', 'W', 'Y', 'F') or
            seq[i+5] in ('L', 'W', 'Y', 'F')
    ):
        print("%s\n%s\n%s" %
            ("".join(seq[0:start]),
             "".join(seq[start:start+7]),
             "".join(seq[start+7:])
            )
        )

This will produce the following output when run:

$ foo.py
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMV
PKCFNEL
LMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

If you insist on using regular expressions, you could separate it into two expressions and use something like this:

import re
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
re1='P[^RKH][^P][^R][LWYF].[^P]'
re2='P[^P][^R][LWYF].[^P]'
r=re.findall('(?:P[^RKH][^P][^R][LWYF].[^P])|(?:P[^P][^R][LWYF].[^P])',X)
print(r)

Before you do, however, please consider the wise words of XKCD 208:


Is there a tool to convert a variableStep .wig to fixedStep one?

This is a temporary solution I came up with that worked for my specific case:

def convert_line(line, chrom, span, start=None, fixedStep=None, count=None):
        new_line=''
        if fixedStep is None or count is None or start is None:
                line_arr  = line.split()
                start = line_arr[0]
                end   = int(line_arr[0]) + int(span)
                val   = line_arr[1]
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)
        else:
                line_arr  = line.split()
                val   = np.float(line_arr[0])
                start = int(start) + int(fixedStep)*int(count)
                end   = start + int(fixedStep)
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)

        return new_line

def convert_file(input_file, output_file):
        chrom=None
        span=None
        fixedStep=None
        count=None
        start=None
        with open(input_file, 'r') as fin, open(output_file, 'w') as fout:
                for line in fin:
                        if len(line)<1 or line[0] == '#':
                                continue
                        elif 'variableStep' in line or 'fixedStep' in line:
                                for key in line.split():
                                        if 'chrom=' in key:
                                                chrom = key.split('chrom=')[-1]
                                        if 'span=' in key:
                                                span = key.split('span=')[-1]
                                        if 'step=' in key:
                                                fixedStep = key.split('step=')[-1]
                                        if 'start=' in key:
                                                start = key.split('start=')[-1]
                                                count = 0
                        else:
                                line_out = convert_line(line, chrom, span, start, fixedStep, count)
                                fout.write(line_out)
                                if count is not None:
                                        count += 1

input:

variableStep chrom=chr16 span=20
59965   7
59985   10
60005   10

output:

chr16   59965   59985   7
chr16   59985   60005   10
chr16   60005   60025   10
chr16   60025   60045   10


In your dataset2 = list(), the MAF argument should be MAF=2_subset[,12] instead of MAF=1_subset[,12]. Since your two datasets have different numbers of lines, this is what causes the error.

I am not sure whether this is a typo in your code, or it is because you do not have MAF information in the second dataset. What I like to do before COLOC is to filter two datasets so they have the same number of lines, and then feed those filtered datasets to the coloc function.

(Also just out of curiosity, is 1_subset a data frame for your GWAS results? I think it is not recommended to have variable names starting with a number.)


How many types of DNA methylation are there?

If I'm understanding you correctly, by "types of DNA methylation" you mean "nucleotide contexts where DNA methylation occurs".
This is going to be a function of the methyltransferase proteins involved in the methylation process, and this largely depends on the organism, and even cell type being studied.

Generally, what has been observed is:

The phenomenon of DNA methylation (methyl group addition to A/C nucleotides) has been observed in both eukaryotes and prokaryotes
The nucleotide context of DNA methylation differs between the organism, and even the cell type being studied, due to the methyltransferase proteins involved in the methylation process
Methylation of RNA, predominantly in N6 position As

In humans:

The DNMT1 methyltransferase is referred to as the "maintenance" methyltransferase, due to its observation of binding to hemi-methylated CpG sites (methylated C on one strand and not the other)
The DNMT3A/B methyltransferases are referred to as the "de novo" methyltransferases for their affinity of binding to naked DNA
These proteins almost exclusively methylate CpG dinucleotides in human cells
There are cases of non-negligible methylation of cytosines in CHG and CHH nucleotide contexts (brain tissues, for example)

The proteins and their effects on DNA methylation are similar for mice, as well.
Plants tend to have similar nucleotide contexts of methylation, but have much more CHG and CHH methylation than mammals.

In prokaryotes, again, the nucleotide context is dependent on the proteins.
Like you mentioned, we have that:

The Dam methyltransferase methylates the N6 position of the A in GATC
The Dcm methyltransferase methylates the C5 position of the second C in CCAGG and CCTGG
The EcoKI methyltransferase methylates As in AACNNNNNNGTGC and GCACNNNNNNGTT

NEB has a good list of its restriciton enzymes, and how DNA methylation affects their binding/cutting, which also gives some info on the type of methylation in cells.

As for a "comprehensive list", you'll likely have to find some relevant review papers in your organism(s) of interest and compile on yourself.

But also think about how you're phrasing your question.

  Is there a comprehensive list of which sequences of base pairs allow methylations to occur?

I don't think that many papers studying DNA methylation think of the nucleotide context allowing methylation to occur so much as observing the context in which methylation occurs and how it relates to proteins in those organisms.


This is not a very simple problem. I saved the table and "test.txt" and visualize it below:
library(pheatmap)
x = read.delim("test.txt")

head(x)
  Gene.Symbol            Composite_chr_coords  Proneural      Neural
1      CDKN1B   [36.1:12:12761576-12766570:+] 0.03353670  0.03353670
2       EPB41    [36.1:1:29086215-29319545:+] 0.03882308  0.03882308
3        CLGN  [36.1:4:141529057-141568232:-] 2.12796056  0.41999847
4      PDE10A  [36.1:6:165664525-165995575:-] 1.66012916  0.25349464
5     RALGPS2  [36.1:1:176960923-177152465:+] 0.78971500 -0.01945154
6        TAF5 [36.1:10:105117714-105138812:+] 1.08354870  0.09461930
    Classical Mesenchymal
1  0.03353670 -0.84379400
2 -0.84821990  0.03882308
3  0.41999847  0.41999847
4  0.25349464  0.25349464
5 -0.01945154 -0.01945154
6  0.09461930  0.09461930

You can see there are no so called markers genes and if you look at the scale of the data, it was normalized in microarray and it will be very hard for you to get your RNA-seq to fall into that range, and apply the medoids to identify your clusters or type.
Since you did not provide any data, and it is not trivial to simulate data for this, I can suggest the approach in words here:

Subset your count table to these exact set of 840 genes.

Normalize the data and do a PCA or hierarchical clustering to see if you indeed observe 4 or 3 or 2 subtypes

If you indeed see 4 or 3 clusters this would be straight forward, take the centres of means of these clusters, and correlate with each of the sub types here. Chances are the subtype that is most correlated with one cluster, is the subtype you are interested in.


My understanding is that these are cDNA sequences from the RIKEN project that don't yet have another gene name. The project carried out a comprehensive sequencing of the mouse genome, and as a result found a lot of cDNA transcripts that were previously unknown. Now they have some annotation (which can be varying degrees of functional annotation), and an accepted official identifier to allow people to talk about these transcripts when discovered in other studies.

It'd be interesting to see if any novel transcripts have been found from the Tabula Muris project. They would also need their own identifier, possibly leading to similar intrigue about their history.


Sleuth: transcripts with beta close to 0 are considered differentially expressed in a likelihood-ratio test


Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header?

According to PacBio Sequel BAM format manual, you can find information on the chemistry version in the header entries @RG (read group), tag DS (description), entry BASECALLERVERSION .

In case of future very PacBio-specific questions, you can contact Liz Tseng on Twitter (@magdoll), one of the most responsive PacBio developers. They also have a dedicated SeqAnswers forum, which is directly monitored by them and you are very likely to get a fast and accurate answer.


Is it possible to use SNP heterozygosity as a proxy for Indel heterozygosity?

SNPs, insertions, and deletions are generated by different biological processes (see for example here). I would expect the per-generation rates of these mutational processes to vary across species. I would also expect the degree to which each mutational process varies across species to vary as well. I think the answer to your question is "probably not without further information".

However, priors are a feature of Bayesian statistics. In Bayesian inference, it is good practice to redo one's analysis using several priors to establish how sensitive one's results are to the chosen prior. I would suggest that you try to call indels using a range of priors and compare the quality of your call sets to determine how important the choice of prior actually is.


What are the different kinds of bioluminescent genes?

There are three processes wherein light is emitted.

Bioluminescence: a chemical reaction releases light. The enzyme that does this is luciferase, while its substrate is luciferin, a small molecule. There are half a dozen organisms that have analogous (evolved independently) enzymes and totally different substrates, such as fireflies, Alivibrio bacteria, diatomes, deep sea fish, some mushrooms etc. The wiki article has a list of these and their small molecules.
Fluorescence: light is absorbed by the molecule and is emitted at a lower energy (wider wavelength, redshifted). Green fluorescent protein of Aequorea victoria jellyfish is the classic example of this. It is excited with UV and emits green light. Here is an interactive model, I threw together. Variants exists for biotech reasons that span the whole of the visible spectrum. Other species have a homologous beta barrel, such as corals (red).
Phosphorescence: light is absorbed by the molecule and emitted later.  Ignoring the picosecond delay in fluorescent small molecules (cf. fluorescent lifetime microscopy), this is exclusively a properties of lanthanide elements (e.g. Europium). Lathanide binding domains have been designed, but do no exist in nature.


Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern

First option

My solution uses the following steps:

use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C )
use jjs (java scripting engine) and the htsjdk library to build a bufferof reads having the same name. If any read starts with a pattern, print the whole buffer. The htsjdk should be already contained in the picard.jar.
sort back to coordinate order using samtools.

All in one:

java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
jjs -cp /path/to/picard.jar  script.js |\
samtools sort -T TMP -o out.bam -

with script.js

var out=java.lang.System.out;
var File = Java.type("java.io.File");
var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
var iter=samReader.iterator();
var buffer=[];
for(;;)
        {
        var rec=null;
        if(iter.hasNext()) rec=iter.next();
        if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
            {
            var i=0;
            for(i=0;i< buffer.length;++i)
                {
                if(buffer[i].getReadString().match(/^AAAA/)) break;
                }
            if(i!=buffer.length)
                {
                for(i=0;i< buffer.length ;++i)
                    {
                    samWriter.addAlignment(buffer[i]);
                    }
                }
            if(rec==null) break;
            buffer=[];
            }
        buffer.push(rec);
        }

samReader.close();
samWriter.close();

A second option is to use samtools view and awk to export the read names:

 samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt

and to use Picard FilterSamReads and READ_LIST_FILE to filter the reads.

UPDATE: 3rd solution:

I've updated my tool samjdk: there is now a new option --pair to work with query-name-sorted BAM files.

The command would be:

java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
java -jar dist/samjdk.jar --samoutputformat BAM --pair \
     -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
samtools sort -T TMP -o out.bam -


What is the typical host-to-bug DNA ratio found in nanopore sequencing without amplification?

My recollection of a conversation I had with Justin O'Grady is that it's about 99% host DNA for unfiltered samples. There are a few ways to reduce that portion prior to sequencing; he has written a paper that mentions a saponin-based host depletion method that is quite effective:
https://doi.org/10.1038/s41587-019-0156-5


Melt p-values for CpG sites mapping to the same gene

Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.


If using DESeq2 with GSEA, I'd recommend ranking by shrunk log2FC values.

It'd also be worth considering ranking positive and negative associations separately, because the standard GSEA algorithm doesn't cross at the zero point when associations change from positive to negative.

You shouldn't be using p-values to rank anything. A p-value gives information about the reliability of a result, but not its importance. See point #5 of this article about p-values:

https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i31


How can I obtain the percentage gene expression per identity class in Seurat as further processible numbers (e.g. matrix)?

This can be solved like this:

library(Seurat)

my_genes <- c("gene1", "gene2", "gene3")

exp <- FetchData(object, my_genes)

matrix <- as.matrix(colMeans(exp  > 0))*100


How to query NCBI (Nucleotide database) by a feature qualifier?

For this particular query, try host=Bemisia tabaci B[WORD].  However, you may also be interested in related queries like host=Bemisia tabaci B biotype[WORD].
But, I am not sure a general solution exists, as I do not believe all qualifiers are indexed.

Also, you can see how your query is interpreted on the right side under "Search details".  For the original attempt ([host]"Bemisia tabaci B"), it is interpreted as:
host[All Fields] AND "Bemisia tabaci B"[All Fields]

while the goal is to produce
host=Bemisia tabaci B biotype[WORD]


Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option:

  -s, --samples [^]LIST

  Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools view is the exception where some
  tags will be updated (unless the -I, --no-update option is used; see bcftools view documentation). To use updated tags for the subset in another command one can pipe from view into that command. For example:

  -S, --samples-file FILE

  File of sample names to include or exclude if prefixed with "^". One sample per line. See also the note above for the -s, --samples option. The command bcftools call accepts an optional second column indicating ploidy (0, 1 or 2) or sex (as defined by --ploidy, for example "F" or "M"), and can parse also PED files. If the second column is not present, the sex "F" is assumed. With bcftools call -C trio, PED file is expected. File formats examples:

sample1    1
sample2    2
sample3    2

  or

sample1    M
sample2    F
sample3    F

  or a .ped file (here is shown a minimum working example, the first column is ignored and the last indicates sex: 1=male, 2=female):

ignored daughterA fatherA motherA 2
ignored sonB fatherB motherB 1

Example usage:

bcftools view -s sample1,sample2 file.vcf > filtered.vcf
bcftools view -S sample_file.txt file.vcf > filtered.vcf

See the bcftools manpage for more information.


Demultiplex nanopore reads with custom barcodes

My semi-automatic workflow for demultiplexing nanopore reads using LAST can be found here - pay attention to the comments, where I've identified bits in the steps which might trip up people trying to follow them:

https://dx.doi.org/10.17504/protocols.io.7vmhn46

While this protocol is designed around using the ONT barcodes, the barcode fasta file can be replaced by another file [using barcode ID in the header] without changing the protocol. Make sure that this file only includes the unique parts of the barcodes (with a few bases of context if necessary to exclude random matches), so that reads won't be matched to the non-unique part of the wrong barcode sequence.


No, there is currently no GFF support in biopython.

However, you can read in GFF files into python using this package, gffutils. There are also a few other packages to read/write GFF files, like gff3.


Tool for calibrating base quality scores?

I think what you are looking for is GATK's Base Quality Score Recalibration.  Their documentation is great so I won't rewrite it here.

  BQSR stands for Base Quality Score Recalibration. In a nutshell, it is
  a data pre-processing step that detects systematic errors made by the
  sequencing machine when it estimates the accuracy of each base call.


Why are my Chi-squared test results different from those in a published table?

Your calculations seem right, perhaps there was an error on their side.

I also looked at the number of samples reported, but they use the same amount of samples in each case. Because they are testing whether Linc-ZNF469-3 is different in several tissues, they should have used a multiple test correction for the p-value (although reporting the original values is usually also ok).

To answer your question in the comments. You can't correlate a categorical variable with a continuous  variable. Or they did another test or they are not explaining clearly what they did.


Difference between genometools installed with anaconda vs apt-get?

These appear to be two separate collections of tools. The conda version is a collection of python scripts whereas the apt version is based on a C library. So, as you might expect, the python version does not have a gt binary.

They also have separate GitHub accounts conda genometools vs apt genometools.


How does one construct a cladogram of intraspecies relationships?

Even with inbreeding and other genetic phenomena that might mask actual evolution of these cultivars, any phylogenetic methodology would be capable of determining relationships accurately.

Try creating a Neighbour Joining tree with MEGA, which is one of the simplest methods available. This should give you enough to check the relationships of the cultivars.


Calculate molecular mass from chemical formula C9H15N4O8P with Python?

As you mentioned, rdKit does not accept a chemical formula as input. The main reason for this are probably isomers. You could however use rdKit to initiate all the atoms, get their mass and add them together.

import re
from rdkit import Chem

def get_mass(formula):

    parts = re.findall("[A-Z][a-z]?|[0-9]+", formula)
    mass = 0

    for index in range(len(parts)):
        if parts[index].isnumeric():
            continue

        atom = Chem.Atom(parts[index])
        multiplier = int(parts[index + 1]) if len(parts) > index + 1 and parts[index + 1].isnumeric() else 1
        mass += atom.GetMass() * multiplier

    return mass

Then...

>>> get_mass("C9H15N4O8P")
338.21299999999997


How to incorportate RIN values as covariate in the design matrix?

In the DESeq2 manual there is a section titled "How can I include a continuous covariate in the design formula?" that deals with your question.
Basically the process is no different from using a discrete covariate.

DESeq <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ RIN + condition)

In this example condition is the experimental groups and RIN is the continuous covariate for the RNA-integrity.

For this to work properly there has to be a linear relationship between the RNA integrity and the gene expression values. As stated in the manual: "Continuous covariates might make sense in certain experiments, where a constant fold change might be expected for each unit of the covariate."
Therefore, it may be better to code the RIN variable as a factor like RIN_high or RIN_low. This is more robust than using the absolute values of the variable itself.

However, from your desription I get the impession that the RIN values are generally lower in the treatment group compared to the control group. If this is the case, using RIN as a covariate will by of limited use. DESeq2 will not be able to discriminate if a expression difference is due to a difference in RNA integrity or due to your treatment. Also, in this case coding RIN into high and low would result in an error because the RIN and treatment variable would define the same groups.

Generally, covariates work best when they are evenly distributed across the experimental groups.


How to map selected genes to Metabolic pathway Maps

KEGG has a webtool for this, KEGG Mapper.

You can just copy-paste your gene ids in there and see how your pathways are colored.


Is there any tool for fuzzy sequence matching?

Python solution:
e<=2 will allow up to 2 insertions, 2 deletions, or 2 SNPs in any 2 combinations

import regex
seq_list=["GGGAGGG", "GGGAAGGG", "GGGAAAGGG", "TTTATTT", "TTTAATTT", "TTTAAATTT", "GGGG", "GGG"]
for sequence in seq_list:
    if regex.search("(GGGGGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAGGG
GGGAAGGG
TTTATTT
TTTAATTT
GGGG

Can also include a wildcard

seq_list=["GGGAAAAAGGG", "GGGAAAAATGG", "GGGAAAAATTG", "GGGAAAAATTT"]
for sequence in seq_list:
    if regex.search("(GGG.*GGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAAAAAGGG
GGGAAAAATGG
GGGAAAAATTG


To end up with a single table made by merging all your tables, you can give a try with reduce() of purrr:  The first argument .x would be the list of tables to merge (in your case the annotation file plus all the other files) and the .y argument would be your function, i.e. inner_join(). reduce() would merge all the tables within the given list, left to right or right to left. ... would enable you to specify column(s) to use as "anchors" while merging.

If you would prefer to keep your tables separately then lapply() is the answer: lapply(list_of_files, function(x) inner_join(annotation_file, x)). This would give a list of merged tables.


I'm working on something similar as a side project. I'll detail the steps that I'm using, but I'm curious if others can suggest a better way.

You don't mention if your data are paired-end, but in my case I have paired-end sequences that overlap in the middle. I used FLASh to merge the reads and seqtk to convert the reads to fasta format.
I used the cluster_fast algorithm in vsearch to cluster sequences at 100% identity, then filtered out rare sequence types (<10 sequences) to account for sequencing errors.
This resulted in a manageable number of unique sequences, which I searched against a custom database using vsearch, but which could also have been searched against the nr database using blast.


Have a look at BCFtools, which I think does most of the stuff that VCFtools does. I notice there's a --mask option in BCFtools consensus, which may do what you want:
-m, --mask FILE
 BED file or TAB file with regions to be replaced with N
 (the default) or as specified by the next --mask-with
 option. See discussion of --regions-file in Common
 Options for file format details


Can I use BioEdit to quantify the sequence similarity?

No Bioedit is an alignment editor only. It is not in anyway an analysis program.

Dump the format in fasta or mega (its the option is available) and import it into Mega X here. You can quantify the diversity as a pairwise matrix, population analysis (if you have a clear a priori that you can subdivide your sequences into different groups) or phylogeny ... which provides diversity as mutations per alignment position. Its pointy and clicky and easy to use.


Article GEO gives counts that are not integers. Should UMI counts be integers?

The data processing section doesn't mention using UMIs at all, so it's likely there were none (whether they'll be present depends on how the libraries are made, with lower throughput methods being less likely to include them). The values are stated as being expected counts from RSEM, which is why they're not always integers.


PCA plot shows big difference but not many differentially expressed genes are found

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.

On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.


Biological datasets for polythetic classification

I may not be understanding exactly what you are talking about, but it seems to me that nearly any set of classes of interest must be polythetic. Most groups of comparable things differ in some characteristics but also have characteristics in common, and none of those characters is 100% identical with membership in the class.
Some examples off the top of my head:

Biological species / taxonomic groups (classify according to phenotypic characters)
Biological species / taxonomic groups (classify according to molecular sequence e.g. DNA)
Cell types (classify according to cell epitopes)
Cell types (classify according to transcriptional profile)
Ecosystems (classify according to species profile)
Disease state (classify according to any one of many molecular profiles)
Members of population groups (classify according to genetic variation)

etc.


Gene expression Table to Expression Matrix converstion

Your question stated differently is as follows:

I have data in long format and need to change it to wide format. For that, I want to keep the gene_raw, use patient_id for the column names and gene_code for the row names. How can I do this easily?

Assuming you know how to read that into R as a data.frame that I'll call d:
library(tidyr)
library(dplyr)
reformatted = d %>% select(patient_id, gene_code, gene_raw) %>% spread(patient_id, gene_raw) %>% column_to_rownames(gene_code)

Or something along those lines.


Can I change the name of file features.tsv to genes.tsv

You can follow the 10X instructions to load your matrix into R:
the following is a direct quote from their webpage:

library(Matrix)
matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

You can then create a seurat object using mat

object <- CreateSeuratObject(mat, project="")


Cellranger aggr and Seurat merge difference?

Cell Ranger aggregate subsamples reads (unless you select none), so you will end up with less total reads in samples that have more initially. The output is still raw counts, but you will have more or less per cell.

Seurat just merges the raw counts matrices and normalizes those.


Is it wise to use RepeatMasker on prokaryotes?

If I understood correctly your question, you want to mask those regions in a (FASTA?) genome. I think you could identify those regions using mummer and mask them using bedtools.
# align genome against itself
nucmer --maxmatch --nosimplify genome.fasta genome.fasta

# select repeats and convert the corrdinates to bed format
show-coords -r -T -H out.delta | awk '{if ($1 != $3 && $2 != $4) print $0}' | awk '{print $8"\t"$1"\t"$2}' > repeats.bed

# mask those bases with bedtools
bedtools maskfasta -fi genome.fasta -bed repeats.bed -fo masked.fasta

Have a look at nucmer and bedtools maskfasta options to fine-tune your analysis.


I am the developer of Uberon and I would be happy to help you with what you need, either from Uberon, or from the FMA.

You mentioned you need something simpler than FMA. There are a variety of tools for creating custom subsets. Additionally, some ontologies provide ready-made subsets for particular purposes.

For example, we make a subset called 'basic.obo' from Uberon:
http://uberon.github.io/downloads.html#subsets

You can examine this directly in either oboedit or Protege5 (note Protege5 parses obo so you don't need to convert).

From your question, you didn't indicate you needed programmatic access, but if you do there are a variety of options, as well as more programmer-friendly JSON exports.

Most ontologies have a number of different options for web-browsing, you can see some of the options for Uberon here:
http://obofoundry.org/ontology/uberon.html

You mentioned annotating disease sites. There may be a tools (either desktop or web-based) for helping with this - e.g. providing autocomplete over a desired subset of an ontology.

You may also want to ask your question on a mailing list like obo-discuss. Unfortunately many ontology developers don't read stackexchange, but most are always keen to hear of community requirements, to help them provide better and more useful products.


Relaxed and sequential Phylip format conversion

I have nutted out a solution for you.
If this is the content of input.fa:
>virusA_detailed_description
ACGTGACGATCCCCAAACGTGACGATCCCCAA
>virusB_detailed_description
ACGTGACGATCGGGAAACGTGACGATCGGGAA
>virusC_detailed_description
ACGTGACGATCAAAAAACGTGACGATCAAAAA
>virusD_detailed_description
ACGTGACGATTTTTAAACGTGACGATTTTTAA

And this is the script:
#!/usr/bin/env python3

'''
Converts a fasta file to phylip format
'''

from Bio import AlignIO
from Bio.AlignIO.PhylipIO import PhylipWriter
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
idwidth = int(sys.argv[3])
with open(infile, 'r') as input_handle:
    align = AlignIO.read(input_handle, "fasta")
    with open(outfile, 'w') as output_handle:
        PhylipWriter(output_handle).write_alignment(align, id_width=idwidth)

Then, save this script as phylpwrtr.py and do:
chmod +x phylpwrtr.py
./phylpwrtr.py input.fa output_longids.phy 30 # 30 is the width of the header, adjust the value to adjust header length

The contents of output_longids.phy looks like:
 4 32
virusA_detailed_description    ACGTGACGAT CCCCAAACGT GACGATCCCC AA
virusB_detailed_description    ACGTGACGAT CGGGAAACGT GACGATCGGG AA
virusC_detailed_description    ACGTGACGAT CAAAAAACGT GACGATCAAA AA
virusD_detailed_description    ACGTGACGAT TTTTAAACGT GACGATTTTT AA


Imputing missing genotypes from separate genotyping panels

Given that you mention wanting to use 1000 Genomes as a reference panel for imputing genotypes into your two SNP chip panels, I am going to assume that you are working with human data.

In that case there are several options you can go with:

If your two panels are of European descent, then you are probably best off using the HRC reference panel together with a fast genotype imputation tool such as Beagle 4.1 to impute genotypes in each of your two SNP chip panels separately.
If your panels are not of European descent, then you will likely want to use the 1000 Genomes phase 3 reference panel with Beagle 4.1, Impute2, or Minimac3.

In either case, there are two phasing services available that will do much of the heavy lifting for you 1,2.

The second Wellcome Trust Case-Control Consortium paper performed a cross-imputation analysis as you describe.  I don't see many studies using multiple SNP chip panels.  You will need to take care in your analysis that you are not hit by batch effects from using two different SNP chip panels.

Also, none of these methods will work if the region you are imputing into has too few variants.  I'm not sure what the minimum number of variants is, but if you are using a whole genome genotyping panel of at least 500k SNPs, then you should be ok if you impute a whole chromosome at a time.


How can I plot multiple sequence alignment in R as a heatmap

@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.

@StupidWolf:

We read in the alignment:

library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300)

Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :

aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:

bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:

library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.


Why does the Genbank phylogenetic tree of 4000 sars-cov-2 sequences only display 200 nodes?

Okay, I'll provide a formal answer, Genbank phylogeny is to get a gist of the diversity. It isn't a formal phylogenetic analysis. You can down the tree as a pdf, I don't think you can download the treefile (in fact I'm pretty certain you can't). Genbank will perform a reasonable distance method and from recollection it will "collapse clades", so if you see little triangles on the tree it represents a large number of taxa represented as a triangle.

I was surprised to learn there are 4000 SARS-CoV-2 sequences.

To make a tree you need to download the data, align (muscle) and send through a formal phylogeny program. If you think Genbank is okay, then try MegaX its very point and clicky. Generally, maximum likelihood phylogeny has returned to vogue for at least 1 decade now (previously Bayesian). However, the Beast package is a Bayesian analysis which is still going strong and the main tool now used in molecular dating.

In summary, to perform a tree analysis beyond Genbank

use the download option (there are two variants - the entire genome or the cut locus based on the query sequence)
align using muscle or clustalo (muscle is really quick(
perform a tree analysis using MEGA X, phylml, raxml or Beast depending on what you objectives are


What is "block-compressed" file in samtools?

A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block so that the whole file does not need to be decompressed and therefore can be accessed randomly. This isn't true of files compressed as a single block using gzip.

As far as I was aware samtools only works with uncompressed FASTA files.


Finding 2000 bp sequences that occur more than once in the genome

MUMmer includes a repeat-match tool which can do exactly that:
$ repeat-match
USAGE:  repeat-match  [options]  <genome-file>

Find all maximal exact matches in <genome-file>

Options:
 -E    Use exhaustive (slow) search to find matches
 -f    Forward strand only, don't use reverse complement
 -n #  Set minimum exact match length to #
 -t    Only output tandem repeats
 -V #  Set level of verbose (debugging) printing to #

Note that a "maximal exact match" is just a maximal repeat.


R not being able to find the file despite being in correct directory

The problem is just what it says.  By omitting the quotes around "Biobanque.csv", R thinks you are talking about an object whose name is Biobanque.csv. And R is confused because of course you have not created an object by that name.


How can I combine a kinship matrix with subset individuals when using rvtests?

Adam Locke (a collaborator of this project) suggests that removing covariate information for the unselected individuals (i.e. setting it to NA) works around this problem:

  I believe the problem is that he is using a pre-computed kinship matrix including both males and females, and when using the “—peopleIncludeFile” it can’t properly select the right people from the kinship matrix. I don’t know that others have tried to do analysis this way. An alternative, that I think has worked for others, is to generate phenotype variables that have values only for males or females (so, e.g., BMI, BMI_females, BMI_males). For whatever reason, this appears to correctly only use the and select the samples desired without the kinship matrix.

Presumably a kinship matrix that is generated for just the subset individuals will work as well, although I don't like the idea that rvtests isn't able to do a label lookup on the kinship matrix to find the correct individuals.


Not easy to reproduce your error but here is my guess. This is a string, you can check it:
input$dataset
class(input$dataset)

This is a simple solution. If you want to call out the dataset depending on input, use a list:
set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")
data_list = list(set_1=set_1,set_2=set_2)

Then in the server :
datasetInput <- reactive({
        df <- data_list[[input$dataset]]
            })


Ensembl now have a tool called Haplosaurus (still in beta) which will convert your phased VCF into actual protein sequences.

  Haplosaurus

  haplo is a local tool implementation of the same functionality that
  powers the Ensembl transcript haplotypes view. It takes phased
  genotypes from a VCF and constructs a pair of haplotype sequences for
  each overlapped transcript; these sequences are also translated into
  predicted protein haplotype sequences. Each variant haplotype sequence
  is aligned and compared to the reference, and an HGVS-like name is
  constructed representing its differences to the reference.

  This approach offers an advantage over VEP's analysis, which treats
  each input variant independently. By considering the combined change
  contributed by all the variant alleles across a transcript, the
  compound effects the variants may have are correctly accounted for.

  haplo shares much of the same command line functionality with vep, and
  can use VEP caches, Ensembl databases, GFF and GTF files as sources of
  transcript data; all vep command line flags relating to this
  functionality work the same with haplo.

  Usage

  Input data must be a VCF containing phased genotype data for at
  least one individual; no other formats are currently supported.

  When using a VEP cache as the source of transcript annotation, the
  first time you run haplo with a particular cache it will spend some
  time scanning transcript locations in the cache.

./haplo -i input.vcf -o out.txt -cache


Reduce number of transcripts in a highly variable de novo transcriptome assembly

As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

Map RNA-seq reads to Trinity-generated transcripts using Salmon
Use the expression of BUSCO genes to set a credible signal cutoff
Subset transcripts based on this threshold [for us it was 50 counts]
Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:

cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

More details here:

https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9


How can I build a protein network pathway from a gene expression quantification file?

There are multiple ways to do this, and multiple protein interaction databases besides the ones you mentioned, such as BioGRID or IntAct. Interaction databases are different in how interactions are defined, sometimes it can be experimental evidence of interaction, sometimes coexpression, orthology-based predictions, etc.

There is no single solution to your problem. For String-DB you can use their R package STRINGdb.

# string-db
library(STRINGdb)

# retrieve full graph for H. sapiens
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )

# define genes of interest
genes.of.interest <- ...

# get their neighbors
neighbors <- string_db$get_neighbors(genes.of.interest)

# get the subgraph of interest (your genes + their neighbors)
my.subgraph <- string_db$get_subnetwork(c(genes.of.interest, neighbors))

# look how many genes and interactions you have
my.subgraph

If you prefer GUI, you might consider using Cytoscape. They have multiple addons for interaction databases, e.g. for STRINGdb. There you can just provide your gene list and get a network in a few clicks.


Calculate the occurrence of motif in sequences (including overlaps) using Python

This is a good task for Biopython:

from Bio import SeqIO

with open('input_file') as f:
    motif = next(f).rstrip().lower()
    for r in SeqIO.parse(f, 'fasta'):
        print(r.id, ':', r.seq.lower().count_overlap(motif))

Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.


How to detect a mutation and predict its consequence?

I'll follow up to the great answer from Kamil S Jaron:

Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.

If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.


truncate fastq.gz file without error in Snakemake

Snakemake is using bash strict mode to ensure best practice error reporting in shell commands. Unfortunately some tools violate bash strict mode.

More details of the error, here. Alternatively, you can use if you don't want disable the strict mode.

zcat file.fastq.gz  | awk '(NR<=400000)' > subset.fastq


Is there an R package that computes homoplasy excess ratios (HER)?

Maynard's Smith's homoplasy is what you want and it is available in the recombination detection package (RDP) described here. Its really a Windows package rather than R.


RNA-seq gene/transcript read counts database for Mouse

EBI has an expression atlas, including mouse: https://www.ebi.ac.uk/gxa/experiments?species=mus+musculus

If you click on one of the individual experiments you can download tables of TPM or RPKM per gene.

There are currently about 1000 individual experiments in different tissues etc.


Using the blast command line tool, you can limit your search to a specific taxon using the –taxids option.

So, in your case -taxids 3193 should limit your search to plants.

If you are using the ncbi blast webinterface, you can restrict your search by entering "plants" into the Organism field in the "Choose Search Set" section


Using BLAT command line tool to blat split sequences

The FAQ offers an answer:

  I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.

  We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. The web-based blat is tuned permissively with a minimum cut-off score of 20, which will display most of the alignments. Other than to confirm that your command-line blat is working, there is little use in perfectly replicating the web-based blat results. We advise deciding which filtering parameters make the most sense for the experiment or analysis. Often these settings will be different and more stringent than those of the web-based blat. With that in mind, use the following settings to replicate the search results of the web-based blat:

  [ . . . ]

  standalone blat:

  blat search:

blat -stepSize=5 -repMatch=2253 -minScore=0 -minIdentity=0 database.2bit query.fa output.psl


Outputting crunchable list of HBonds from Pymol

Thanks a lot to all who commented/replied, the author of the scripts provided a response. It was down to a syntax change in PyMol of which I wasn't aware that caused my selection to come up empty.


Preparing PDB file for modeling with Swiss-Model

As explained in the website the data introduced must be:

  Format must be FASTA, Clustal, plain string, or a valid UniProtKB AC

So you need to paste just the protein sequence. Without seeing the data you paste I would guess that you don't use a valid plain string.

As Jonathan Moore explains below a FASTA format is:

>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL

(breaking the sequence in several lines is optional but usually done fore easy reading)


How to convert data in gmt format to dataframe?

Ok. Answering my own question.

I did this way. And I got the output I want.

install.packages("msigdbr")
library(msigdbr)

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

Output:

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

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


What is the --rev_comp_mapping_barcode parameter in the QIIME1 script extract_barcodes.py?

As the help page says:

  --rev_comp_mapping_barcodes
      Reverse complement barcode in mapping before lookup (useful if barcodes in mapping file are reverse complements of golay codes) [default: False]

As the help page of extract_barcode says:

  --rev_comp_bc1
  Reverse complement barcode 1 before writing [default: False]

You seem to be doing two complete different things. extract_barcode.py is designed to

  format fastq sequence and barcode data so they are compatible with split_libraries_fastq.py

While the first one is for:

  demultiplexing of Fastq sequence data where barcodes and sequences are contained in two separate fastq files

So it seems that depending on the machine used or the previous processing you are given with the barcode in another file or together. I would advice to ask the person who gave you this data.


How to add cluster name to metadata in Seurat?

To add the column status can try the following:
B@meta.data$status <- status[row.names(B@meta.data),]$status

To rename the clusters you can try the following:
B@meta.data$seurat_clusters <- B@meta.data$status

However, I think the plotting functions in Seurat do not use the cluster information from B@meta.data, rather they use information from B@active.idents.
So, I guess you have to do the following instead:
B@active.idents <-  B@meta.data$status

This should work, but, I will verify and update this answer later.


BinaryCIF vs MMTF formats, which one to choose?

So binaryCIF is newer and is heavily inspired by MMTF.
NGL, the most "loved" JS library (IMO) to show protein fetches proteins as MMTF, while Mol* uses bCIF —Mol* has the official repo of the bCIF format. Alex Rose, developed the former at the RCSB PDB, and the latter in collaboration with the PDBe. He also worked with Antony Bradley, the main author on the MMTF format, who left the RCSB PDB, some 4+ years ago. So there is definitely no competition there.

The MMTF download of the RCSB works well: https://mmtf.rcsb.org/download.html
The MMTF format did not catch on —okay, mmCIF is only slowly replacing PDB, which is still the only accepted format for analyses—, hence the need for a bCIF format
The libraries to do the conversion aren't without glitches
In the distant future bCIF will be the sole binary format
mmCIF <--> bCIF conversion is easier

However, I really ought to stress that mmCIF is a better choice to bCIF or mmTF. And even then many tools still hate mmCIFs and use PDB format.
If memory footprint is not really a problem, then stick with mmCIF. If you are running a webserver, you may be sending data as gzipped making this a moot point.


Why are my kallisto and salmon results differing so much just for lncRNA transcripts?

There has been a detailed comparison of kallisto and salmon for lncRNA transcripts by Zheng et al https://www.biorxiv.org/content/biorxiv/early/2018/01/09/241869.full.pdf

It shows them to be near identical. This suggests there is some problem with the way you made the comparison/plot.


It looks like rdkit and rdkit's Chem just use the matplotlib. So, if you import matplotlib you can edit your figure like any other matplotlib figure.

I produced this from the command line, so it may work differently in the Jupyter environment.

It looks like Chem doesn't properly scale the drawing based on how large the chemical is, so I added plt.xlim and plt.ylim to make the molecule larger in the frame.

from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
plt.title('Test')
plt.axis('off')
coord_min = 0.35
coord_max = 0.65
plt.xlim((coord_min, coord_max))
plt.ylim((coord_min, coord_max))
fig.savefig('mol.jpeg', bbox_inches='tight')


Orthogonal DNA sequences are the shortest possible single stranded nucleotide stretches constructed. These sequences are used as probes in DNA hybridization. Because of its length and construction it has least chance of binding with other probes.

This technique improves the efficacy and specificity and is clinically used to detect genetic disorders in shortest time duration.

Via:
https://www.nature.com/articles/srep24650

https://www.pnas.org/content/106/7/2289


How can I obtain a list of all NCBI gene ID's along with their full name, symbol, and also known as symbols?

All of these data are in the gene_info.gz file on the NCBI FTP site: https://ftp.ncbi.nlm.nih.gov/gene/DATA/
Specifically, you may want to look at the gene_info.gz file. It has the following fields (shown for an example gene below):
                              #tax_id [  1]: 9606
                               GeneID [  2]: 5768
                               Symbol [  3]: QSOX1
                             LocusTag [  4]: -
                             Synonyms [  5]: Q6|QSCN6
                              dbXrefs [  6]: MIM:603120|HGNC:HGNC:9756|Ensembl:ENSG00000116260
                           chromosome [  7]: 1
                         map_location [  8]: 1q25.2
                          description [  9]: quiescin sulfhydryl oxidase 1
                         type_of_gene [ 10]: protein-coding
   Symbol_from_nomenclature_authority [ 11]: QSOX1
Full_name_from_nomenclature_authority [ 12]: quiescin sulfhydryl oxidase 1
                  Nomenclature_status [ 13]: O
                   Other_designations [ 14]: sulfhydryl oxidase 1|quiescin Q6 sulfhydryl oxidase 1|testis tissue sperm-binding protein Li 62n|thiol oxidase 1
                    Modification_date [ 15]: 20200604
                         Feature_type [ 16]: -

If you are starting  with a set of NCBI GeneIDs and want this sort of information for just those, you can use Entrez Direct for that. If you are starting with RefSeq transcript or protein accessions, be sure to check out the other files in the FTP path mentioned above.


Issue with visualising cladogram/phylogenetic tree with multiple sequence alignment data in R?

You were very close, you just need to supply the offset argument to msaplot().

So using your code block above with one extra argument in msaplot()

library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta", offset = 5)

yields

Thank you for supplying a minimal and reproducible example!


How do PCR duplicates arise and why is it important to remove them for NGS analysis?

In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.

For example, consider the following hypothetical scenario.

                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
               ..........................C.
                ............................
                      .............T..............
                       ..................C.........
                           ............................
                               ............................
                                   ............................
                                    .....C......................
                                       ..C.........................

The marked position has a coverage of 9 (9 overlapping reads), 4 of which suggest the presence of an alternate allele. Since these 4 reads map to different positions, they are independent observations supporting the same putative single nucleotide variant (SNV). The T in the third read can safely be ignored as a sequencing error, since it only appears once.

Now consider the following hypothetical scenario.

                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
                ............................
                      .............T..............
                       ..................C.........
                       ..................C.........
                       ..................C.........
                       ..................C.........
                               ............................
                               ............................
                               ............................

In this case, the marked position also has a coverage of 9, but 7 of those reads appear to be PCR duplicates. If we remove the duplicates, the marked position has only a coverage of 4, and only one C at the marked position, which is not enough confidence to call the SNV.


awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa

An explanation of the non-obvious bits:

BEGIN{FS="::"} split columns using :: as the designator
$1~">" If there's a > in the line
substr($1,2) Trim the > that was at the beginning of the line.


tools to reconcile experimental transcripts with reference annotation

I've never tried this myself, so I don't know how easy this is...

One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You should then be able to get the information you want from the transfrag codes.


How to get Nanopore MinION fast5 from SRA

When I was trying to upload MinION data to SRA, they wouldn't allow me to upload only the raw signal data (i.e. FAST5 files). Every file needed to contain a FASTQ sequence, and this sequence was extracted from the file and stored in the archive. Given this, I don't expect it will be possible to retrieve signal or event information from SRA. I also don't expect this will be possible in the future either, given that the SRA database expects every record to have a FASTQ-like sequence attached.

However, EBI/ENA do allow the upload of raw FAST5 files to the database (e.g. here), but there are no FAST5 reads for the specific single-cell B1a study you have indicated. If you're interested specifically in that study, then I recommend that you ask Miten Jain directly for access to those reads; he should have another method to give you access to them.

ENA is "ready" for ONT reads, but researchers need to upload their reads to ENA in order for that database to be useful.


Frequency of specific viral sequence in .BAM or .fastq

Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.


What would be a reasonable healthy data set/control group to compare with a cell line from TCGA?

In the TCGA datasets there is a variable (type) containing the information if the sample was from the tumor or from an adjacent region, which is usually considered as healthy.

As this samples come from the same patient and the same tissues they are used in the analysis for the comparisons and classifications.


Converting Gene Symbol to Ensembl ID in R

You can use the standard annotation package for humans (you can also use biomaRt, but it can be more confusing, see below):

library("org.Hs.eg.db") # remember to install it if you don't have it already
symbols <- mapIds(org.Hs.eg.db, keys = ensemblsIDS, keytype = "ENSEMBL", column="SYMBOL")

The ensemblsIDS should be a character vector: c("ENSG001", "ENSG002", ..)

The solution with biomaRt:

library("biomaRt")
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
getBM(attributes='hgnc_symbol',
      filters = 'ensembl_gene_id',
      values = ensemblsIDS,
      mart = ensembl)


DESeq2 - Error in lfcShrink: 'coef' should specify same coefficient as in results 'res'


Where can I find .cdf files for Affymetrix ChIP-chip data analysis?

It looks like it's in GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10187

There is apparently a zipped copy linked at the bottom in the supplementary files.

When people upload a study based on an array, there is a corresponding 'platform' record where the array details should be deposited.

Edit: you will probably need to gunzip the file first, to get it into CDF format.


How to convert a Pileup file to VCF format with Hg19 alignment

Short answer

You'll have better luck searching for the term variant calling than with the term converting. A Google search for "variant calling" gives several beginner-level tutorials on the topic. Commonly used tools for variant calling include samtools, GATK, and FreeBayes.

Long answer

Converting data files from one format to another is a common task in bioinformatics. But that is not the task you want to perform here—the term convert is misleading.

A pileup file contains reads aligned to a reference sequence. A VCF file contains a description of putative or validated genetic variant. The process of creating a VCF from alignments is commonly called variant discovery, variant prediction, or variant calling. This is a broad and active area of research interest in bioinformatics, with many available algorithms and tools, each with a different set of parameters and assumptions.

This forum is not well designed to give you a good introduction to this topic. We are much better equipped to answer specific questions you have about how a particular tool or algorithm works, or why you are observing a particular unexpected result, and so on...assuming you have put in some time to research the question yourself.


Voom function from limma package and Normalization on counts data

I'd like to know answer for this myself. It seems the normalize.method argument is set to none by default. You can edit by looking at the methods in normalizeBetweenArrays.
But how do I know if and when I have to use normalize.method in my analysis?
The output from voom is log2CPM of count values (adding 0.5 to them).


Does the MAPQ=0 fraction of a BAM file depend on the insert sizes?

Yes, as a general rule of thumb mappability increases with insert size (up to a limit) and read length. Whether this will actually occur in a given case will depend more on how randomly the sequencing samples from the genome to begin with (i.e., if the library prep happens to select for/against high mappability regions then the insert size won't much matter). The basic reason behind this is that if one end of a read aligns to a repetitive element that the likelihood of being able to use the other end as an anchor will increase a bit if it's further away. Now having said that, Illumina sequencers in particular like a certain fragment length for efficient sequencing, so you can't go too crazy on increasing the insert sizes and still get decent output (unless you modify your library prep).


The UCSC genome browser supports local mirrors. You could make a local mirror that has the "TransMap Alignments" cross-species gene track. From the track description page:

  These tracks contain cDNA and gene alignments produced by the TransMap cross-species alignment algorithm from other vertebrate species in the UCSC Genome Browser. For closer evolutionary distances, the alignments are created using syntenically filtered LASTZ or BLASTZ alignment chains, resulting in a prediction of the orthologous genes in human. For more distant organisms, reciprocal best alignments are used.

They provide a helpful tutorial on setting up a local installation here.


How to construct the phylogenetic tree of the whole genus?

The MLST genes are not used to contruct a tree, because each pattern of alleles determines a predetermined genotype, e.g.  Staphlococcus aureus ST 22 is the MLST type concerning the MRSA epidemic. You could build a tree, but its only purpose is to visualise the ST, so its a bit pointless but okay at least you will not mistype something.

Performing a whole genome tree on Strep for the genus is complicated but doable and there are missing genes and stuff and if you just want typing then its not needed.

So if typing is your goal I would submit your MLST loci to the Steptomyces MLST database,

https://pubmlst.org/bigsdb?db=pubmlst_streptomyces_seqdef


These values come from a microarray so they represent the number of sequences that where attached to the microarray's probes (1007_s_at, ...).

In general yes, a higher value means more expression. The normalizations (GC-RMA) applied are to correct for several effects. They help for testing the probability of the given probe to be a false positive.


Based on your description I think you should have a look at a technique called 'biclustering'.

The example on this page defines the goal of this technique as 'Finding subgroups of rows and columns which are as similar as possible to each other and as different as possible to the rest.'

Since your examples are python-based, you could check out scikit-learn's implementations of biclustering.


Is there a way to save a full pdb file in PYMOL after uploading biological unit?

So on the GUI form of PyMOL on the right hand side you have a list of what is ambiguously called objects in PyMOL (models elsewhere), which can be enabled or disabled by clicking on their name.
The name also has a 1/1 in green or similar. This indicates the state. In your case you have 2 states.

To split a multistate object use

split_state prot_A

To combine two objects one uses the command

create new_name, prot_A or prot_B

Now, this has two caveats which if violated result in a nasty piece of modern art:

if the chains, segment ids and residues are the same
if you had two overlapping macromolecules

To deal with the first, change the chain id on one of them.

alter prot_B, chain='B'
   sort

Segment identifier is a weird thing that most comp-biochem tools dislike, so I like to remove that

alter prot_B, segi=''
   sort

Now when you go to save your new_name object, in the popup of the GUI you get a dropdown menu defaulted to enabled (visible objects that have not be clicked to be disabled). There you can choose new_name.

Aligning two

If someone that is in not the OP author has a template protein and needs to make a dimer of another based on that, the process is the same.

create twin, prot_B
align twin, prot_A and chain B
alter twin, chain='B'
sort

SDF files

If someone else has a multicompound sdf where the state is the header name, I strongly advise not to use split_state. Yes, the python pymol module (not the app) handles states in a really weird way, but states can be selected with state 1 selection algebra and state is an argument for most commands that do not accept state selector.


How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?

Presuming you have deepTools installed:

#/!usr/bin/env python
from deeptoolsintervals.parse import GTF
import pysam

anno = GTF(["some/gtf/or/bed/file", "or/more/than/one/if you like/"])
bam = pysam.AlignmentFile("something.bam")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    e = s + 1
    overlaps = anno.findOverlaps(read.reference_name s, e, matchType=4)
    # For a GTF, parse the results to see if this is the correct type.
    # For a BED3/BED6 file there will be no exons, so anything other than
    # None indicates an appropriate overlap

    do something

deeptoolsintervals in the GTF/BED parsing and custom interval tree library used by deepTools. anno holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the matchType=4 is the pivotal part since that specifies that only overlaps having the exact same start position should be included. Assuming you have a stranded library, you might instead do something like switch between matchType=4 (the same start position) and matchType=5 (the same end position) depending on the orientation of the alignment and the the strand= of the gene.

Note that you'll need to make a set() to hold read names to find the mates (that or jump around the file to get them, which is probably slower).


How do programs like STRUCTURE determine populations from alignments?

structure works with SNP/genotype information, so won't work directly off BAM/SAM files. You would need to convert to VCF (e.g. using mpileup/vcftools), then process using something like PLINK to convert it into a format that STRUCTURE recognises.

Here's the description I have about structure from my PhD thesis. It might be a bit dated now, but I expect most of it is still applicable:

A clustering program called structure was created by Falush et al. (2003) that is used to determine population structure using a set of unlinked genotypes. The expected use case for structure is in the analysis of population structure, but it can also be useful for determining the utility of a given set of markers for categorising (or quantifying) a given trait.

The main benefit of using structure over other methods is that it is able to derive estimated ancestry coefficients from the data, rather than a discrete yes/no for each individual. This allows for the use of cutoff values for the assignment of individuals (useful for tweaking false positive and false negative rates in diagnostic tests), and a quantifiable value where a classification is not useful or does not make sense (as in the case of a continuous trait with high heritability such as height).

The program begins by assigning each individual to a given population. This can be either given in the input data as a ‘popinfo’ flag, or determined randomly by structure. After this, the program constructs a probability distribution for the data, considering the population assignments given, and uses that to estimate population allele frequencies. The program then constructs another probability distribution for these data, and uses that to estimate the population that each individual is likely to have originated from. This process is repeated many times (such a process is known as a random walk), and the underlying theory suggests that the progression of population assignments through the random walk process can be used as an approximation of the actual populations of origin for each individual.

The documentation included with structure describes the general algorithm for each iteration (of maybe hundreds of thousands of iterations) as follows:

Sample (population allele frequencies for iteration m) from the probability of (allele frequency) given (genotypes of sampled individuals) and (population of origin of individuals for iteration m − 1).
Sample (population of origin of individuals for iteration m) from the
probability of (population of origin) given (genotypes of sampled individuals) and (allele frequencies for iteration m).

In more simple terms, the steps are similar to the following statements:

these populations probably have this allele frequency distribution
therefore...
these individuals probably came from this population

The structure program is able to accommodate admixture into its model
by considering the population of origin as a probability, rather than a single
discrete value.

The file that is output from structure is a text file with individual identifiers and predicted ancestry coefficients (Q values). While the program does provide means to display these data graphically, the visual representation of data (including sorting) is better controlled using an external program. The most commonly used program is Noah Rosenberg’s distruct (Rosenberg, 2004), a program that allows colour customisation and group re-ordering. For the purpose of this thesis, a custom R script, snpchip2structure.R, has been designed which also allows selection of different sorting methods, a scatter plot for two populations, error bars, and a few additional features (see figure 1.7). More details about this script can be found in the programs directory of the supplementary CD for this thesis.