Friday 19 February 2016

Extract phenotypes/GO terms for C. elegans genes from WormBase

I had a list of C. elegans genes (WBGene identifiers) and wanted to know which C. elegans genes have known phenotypes. On the WormBase page, they know say that WormMart is retired, and to use WormMine instead. So I did this:

To get phenotypes for a list of C. elegans genes:
1. Went to WormMine
2. Clicked on 'Lists' and uploaded my list of worm genes (separated by commas), and saved the list as 'wormsinglecopygenes2'
3. Then clicked on 'Phenotypes', and chose 'Genes->Phenotypes', and selected the box 'Constrain to be in the list wormsinglecopygenes2', and then pressed 'Go'.
4. Voila! It gave me a nice output spreadsheet.

To get all C. elegans genes with a particular phenotype:
1. Went to WormMine
2. Clicked on 'Templates', and then clicked on 'Phenotype-->Genes'.
3. Set 'phenotype == lethal'.
4. This found 1343 rows. I clicked 'Download'.

Note 12-Aug-2016: I notice that when I use WormMine I get far fewer genes with a particular phenotype than when I search WormBase directly. For example, for the phenotype 'molt defect' ('WBPhenotype:0000638'), when I search using WBPhenotype:0000638 on www.wormbase.org, and download all the results, the results file has 242 distinct genes. In contrast, when I go to WormMine, and choose the template Phenotype->Genes and enter 'molt defect', it gives me a list of just 19 genes. It seems to be missing many genes, eg. WBGene00000039 is found using the website, but not wormmine, and when I checked the wormbase page for that gene, it has an RNAi phenotype observed of 'molt defect' for WBGene00000039, so surely should be listed under WormMine too. I emailed WormBase and heard that RNAi information is not yet used in WormMine however, and that I can get a full file of phenotypes for C. elegans genes from the ftp site: ftp://ftp.wormbase.org/pub/wormbase/releases/WS253/ONTOLOGY/phenotype_association.WS253.wb
(Thanks Michael Paulini for help!)

To get all GO terms for all C. elegans genes:
1. Went to WormMine.
2. Clicked on 'QueryBuilder', then selected 'gene' from the dropdown list.
3. In the browser at the left, for 'Gene'->'WormBase Gene ID', clicked on 'Show'.
4. In the browser at the left, clicked on 'GO Annotations' -> 'Ontology Term' -> 'Ontology Annotations' -> and for 'Identifier' clicked on 'Show'.
5. Clicked 'Show Results' at the bottom right. This gave 274,273 rows of results. Then clicked 'Download'.

Note 7-Jul-2016:
- the website for WormMine now seems to be here.

Thursday 18 February 2016

Calculating distances between drug molecules

Cluster a set of compounds by structural similarity

Given a set of drug molecules from Chembl, I need to calculate pairwise distances between the drugs, for clustering. Noel O'Blog helped me do it using RDKit by Greg Landrum, by writing the Python script below. Given Chembl ids. for the input molecules, it calculates pairwise similarities between drugs based on calculating an ECFP4 fingerprint for each drug molecule, and then calculating the distances between the fingerprints for each pair of drugs. The pairwise similarities given in the output file are scaled from 0-1, so can be converted to distances using 1 - similarity. The input file contains duplicates that are distinguished only by undefined/defined stereocentres; these are likely the same molecule, so the script removes duplicates that are distinguished only by stereochemistry.

Thanks Noel!

import collections

from rdkit import Chem, DataStructs
from rdkit.Chem import Descriptors, AllChem

smiles_lookup = dict( (y, x) for (x, y) in (z.split() for z in open("C:\Tools\Topliss\chembl20\chembl_20.ism") if len(z.split())==2))

def readdata(filename):
    """Return dictionary from chemblid to SMILES"""
    it = open(filename)
    header = next(it)
    data = collections.defaultdict(set)
    for line in it:
        broken = line.rstrip().split()
        if broken[2] in data: # Check for duplicates
            assert False
        data[broken[2]] = broken[1]
    return data

def getmol(smi):
    """Keep a single largest component if disconnected"""
    mol = Chem.MolFromSmiles(smi)
    if "." in smi:
        frags = list(Chem.GetMolFrags(mol, asMols=True))
        frags.sort(key=lambda x:x.GetNumHeavyAtoms(), reverse=True)
        mol = frags[0]
    return mol

# The file contains duplicates that are distinguished only by undefined/defined stereocenters
# These are likely the same molecule. Let's remove duplicates that are distinguished only
# by stereochemistry
def mergeStereoisomers(chemblids):
    smis = [data[chemblid] for chemblid in chemblids_filtered]
    mols = [getmol(smi) for smi in smis]

    normalised = {}
    for i, mol in enumerate(mols):
        smi = Chem.MolToSmiles(mol)
        normalised[smi] = chemblids[i] # Stores only one

    return normalised.values()

if __name__ == "__main__":
    # ECFP setup
    fingerprinter = lambda m: AllChem.GetMorganFingerprintAsBitVect(m, 2, nBits=16384)

    srcs = ["SMILES4CMPDS_smiles.tsv"]
    for src in srcs:
        filename = src
        print "Looking at %s" % filename
        data = readdata(filename)

        chemblids = data.keys()
        print "Distinct chemblids in original TSV", len(chemblids)

        # filter out those molecules without structures in ChEMBL
        # ...this step doesn't apply here
        chemblids_filtered = chemblids

        print "Chembl Ids associated with SMILES", len(chemblids_filtered)
        chemblids_final = mergeStereoisomers(chemblids_filtered)
        print "Chembl Ids after removing stereoisomers", len(chemblids_final)

        smis = [data[chemblid] for chemblid in chemblids_final]
        mols = [getmol(smi) for smi in smis]
        fps = [fingerprinter(mol) for mol in mols]
        N = len(chemblids_final)
        with open("%s_pairwise_sim.txt" % src, "w") as f:
            # header
            f.write("\t".join(chemblids_final) + "\n")
            for i in range(N):
                row = []
                fpA = fps[i]

                for j in range(N):
                    if i==j:
                        sim = 1.0
                    else:
                        fpB = fps[j]
                        sim = DataStructs.FingerprintSimilarity(fpA, fpB)
                    row.append("%.2f" % sim)
                f.write(chemblids_final[i] + "\t" + "\t".join(row) + "\n")



Find any compounds in a new set ('set 2') that are not very similar to those in my original set ('set 1')
A second problem I had was to see if any compounds that I had found from a PubChem search ('set 2') were very different from a list I had already curated from various sources ('set 1'). Noel O'Blog came to the rescue, helping me to use RDKit to find molecules that are structurally different, defined by calculating similarities of <=0.4, when similarities were calculated using the Tanimoto coefficient of ECFP4 fingerprints as implemented in the RDKit toolkit, http://www.rdkit.org). Here's the Python script, which takes 'set 2' in input file 'pubchem_219_compounds', and 'set 1' in input file 'tmp.csv':

import csv

from rdkit import Chem, DataStructs
from rdkit.Chem import Descriptors, AllChem

def readPCCompounds():
    return [x.rstrip().split("\t") for x in open("pubchem_219_compounds")]

def isDigit(char):
    return char >= '0' and char <= '9'

def readX20():
    ans = []
    reader = csv.reader(open("tmp.csv"))

    # Skip header
    for N in range(7):
        reader.next()

    for row in reader:
        if row[0] == "": break
        if row[6][0]=='-': continue
        if row[6].startswith("!"):
            broken = row[6][1:].split("!")
            if len(broken) == 2:
                smiles = broken
            else:
                smiles = [broken[i] for i in range(1, len(broken), 2)]
            for x in smiles:
                ans.append( [row[0], x] )
        else:
            ans.append( [row[0], row[6]] )
    return ans

def getmol(smi):
    """Keep a single largest component if disconnected"""
    mol = Chem.MolFromSmiles(smi)
    if "." in smi:
        frags = list(Chem.GetMolFrags(mol, asMols=True))
        frags.sort(key=lambda x:x.GetNumHeavyAtoms(), reverse=True)
        mol = frags[0]
    return mol

def getFP(smi):
    mol = getmol(smi)
    return AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=16384)

if __name__ == "__main__":
    pubchems = readPCCompounds()
    knowns = readX20()

    pubchemfps = [getFP(smi) for name, smi in pubchems]
    knownfps = [getFP(smi) for name, smi in knowns]

    # for each PubChem entry...
    for pubchem, pubchemfp in zip(pubchems, pubchemfps):
        print "Considering PubChem entry %s..." % pubchem[0],
        # ...compare to each known...
        results = []
        for known, knownfp in zip(knowns, knownfps):
            sim = DataStructs.FingerprintSimilarity(pubchemfp, knownfp)
            if sim > 0.4:
                results.append( (sim, known[0]) )
        if results:
            results.sort(reverse=True)
            print ", ".join("%s (%.2f)" % (x[1], x[0]) for x in results)
        else:
            print "nothing found with sim > 0.4"