\

# Accessing databases

## NCBI

EUtils is a web service offered by the NCBI to access the sequence, literature and other databases by a special format of URLs. PyCogent offers an interface to construct the URLs and retrieve the results in text format.

### From Pubmed

#### Retrieving PubMed records from NCBI by PubMed ID

The process for getting PubMed records by PubMed ID (PMID) is very similar to that for getting sequences. Basically, you just need to pass in the unique id associated with the article. For example, if you want to get the reference to the original PyCogent paper to see how far we’ve come since then, you can do this:

    >>> from cogent.db.ncbi import EFetch
    >>> ef = EFetch(id='17708774', db='pubmed', rettype='brief')
    >>> ef.read()
    '\n1: Knight R et al. PyCogent: a toolkit for makin...[PMID: 17708774] \n'

If you want more information, there are other rettypes, e.g.

    >>> ef = EFetch(id='17708774', db='pubmed', rettype='citation')
    >>> ef.read()
    "\n1: Genome Biol. 2007;8(8):R171. \n\nPyCogent: a toolkit for making sense from sequence.\n\nKnight R, Maxwell P, Birmingham A, Carnes J, Caporaso JG, Easton BC, Eaton M,\nHamady M, Lindsay H, Liu Z, Lozupone C, McDonald D, Robeson M, Sammut R, Smit S,\nWakefield MJ, Widmann J, Wikman S, Wilson S, Ying H, Huttley GA.\n\nDepartment of Chemistry and Biochemistry, University of Colorado, Boulder,\nColorado, USA. rob@spot.colorado.edu\n\nWe have implemented in Python the COmparative GENomic Toolkit, a fully\nintegrated and thoroughly tested framework for novel probabilistic analyses of\nbiological sequences, devising workflows, and generating publication quality\ngraphics. PyCogent includes connectors to remote databases, built-in generalized\nprobabilistic techniques for working with biological sequences, and controllers\nfor third-party applications. The toolkit takes advantage of parallel\narchitectures and runs on a range of hardware and operating systems, and is\navailable under the general public license from\nhttp://sourceforge.net/projects/pycogent.\n\nPublication Types:\n    Evaluation Studies\n    Research Support, N.I.H., Extramural\n    Research Support, Non-U.S. Gov't\n\nMeSH Terms:\n    Animals\n    BRCA1 Protein/genetics\n    Databases, Genetic\n    Genomics/methods*\n    Humans\n    Phylogeny\n    Protein Conformation\n    Proteobacteria/classification\n    Proteobacteria/genetics\n    Sequence Analysis/methods*\n    Software*\n    von Willebrand Factor/chemistry\n    von Willebrand Factor/genetics\n\nSubstances:\n    BRCA1 Protein\n    von Willebrand Factor\n\nPMID: 17708774 [PubMed - indexed for MEDLINE]\n"

Similarly, if you want something more machine-readable (but quite a lot less human-readable), you can specify XML in the retmode:

    >>> ef = EFetch(id='17708774', db='pubmed', rettype='citation', retmode='xml')
    >>> cite = ef.read()
    >>> for line in cite.splitlines()[:5]:
    ...     print line
    ...
    <?xml version="1.0"?>
    <!DOCTYPE PubmedArticleSet PUBLIC "-//NLM//DTD PubMedArticle, 1st January 2011//EN" "http://www.ncbi.nlm.nih.gov/entrez/query/DTD/pubmed_110101.dtd">
    <PubmedArticleSet>
    <PubmedArticle>
        <MedlineCitation Owner="NLM" Status="MEDLINE">

Only a partial example is shown as there are quite a few lines. As with sequences, you can retrieve multiple accessions at once.

#### Searching for records using EUtils

Getting records by their primary identifiers is very nice if you actually have the primary identifiers, but what if you don’t? For example, what if you want to do a search based on a keyword, or have a genbank accession number rather than a gi, or want to get a range of records?

Fortunately, the more general EUtils class allows this kind of complex workflow with relatively little intervention. For example, if you want to search for articles that mention PyCogent:

    >>> from cogent.db.ncbi import EUtils
    >>> eu = EUtils(db='pubmed', rettype='brief')
    >>> res = eu['PyCogent']
    >>> print res.read()

    1: Smit S et al. From knotted to nested RNA st...[PMID: 18230758]

    2: Knight R et al. PyCogent: a toolkit for makin...[PMID: 17708774]

...or perhaps you want only the ones with PyCogent in the title, in which case you can use any qualifier that NCBI supports:

    >>> res = eu['PyCogent[ti]']
    >>> print res.read()

    1: Knight R et al. PyCogent: a toolkit for makin...[PMID: 17708774]

The NCBI-supported list of field qualifiers, and lots of documentation generally on how to do pubmed queries, is here.

One especially useful feature is the ability to get a list of primary identifiers matching a query. You do this by setting `rettype='uilist'` (not idlist any more, so again you may need to update old code examples). For example:

    >>> eu = EUtils(db='pubmed', rettype='uilist')
    >>> res = eu['PyCogent']
    >>> print res.read()
    18230758
    17708774

This is especially useful when you want to do a bunch of queries (whether for journal articles, as shown here, or for sequences), combine the results, then download the actual unique records only once. You could of course do this with an incredibly complex single query, but good luck debugging that query...

### For sequences

#### Fetching FASTA or Genpept sequences from NCBI using EFetch with GI’s

If you already have a list of GI’s (the numeric identifiers that are used by GenBank internally as identifers), your job is very easy: you just need to use `EFetch` to retrieve the corresponding records. In general, this works for any case where the identifiers you have are the primary keys, e.g. for PubMed you use the PubMed ID (see example below).

Here is an example of getting the nucleotide record that corresponds to one particular id, in this case id # 459567 (chosen arbitrarily). The record arrives as a file-like object that can be read; in this case, we are looking at each line and printing the first 40 characters.

    >>> from cogent.db.ncbi import EFetch
    >>> ef = EFetch(id='459567', rettype='fasta')
    >>> lines = ef.read().splitlines()
    >>> for line in lines:
    ...     print line[:40]
    ...
    >gi|459567|dbj|D28543.1|HPCNS5PC Hepatit
    GAGCACGACATCTACCAATGTTGCCAACTGAACCCAGAGG
    GGCTTTACCTTGGTGGTCCCATGTTTAACTCGCGAGGTCA
    CGGGGTTCTTCCAACCAGCATGGGCAATACCCTCACATGT
    GCAGGCCTCACCAATTCTGACATGTTGGTTTGCGGAGATG
    TC

Similarly, if your id refers to a protein record, you can get that by setting the `rettype` to ‘gp’.

    >>> genpept = EFetch(id='1234567,459567', rettype='gp').read()

You’ll probably notice that the lines look suspiciously like FASTA-format records. This is in fact true: the `rettype` parameter controls what type of record you get back. For example, if we do the same search with `rettype='brief'`, we get.

    >>> from cogent.db.ncbi import EFetch
    >>> ef = EFetch(id='459567', rettype='brief')
    >>> lines = ef.read().splitlines()
    >>> for line in lines:
    ...     print line
    ...
    D28543 Hepatitis C virus... [gi:459567]

The current `rettypes` (as of this writing on 4/14/2010) for the ‘core’ NCBI databases are native, fasta, gb, gp, gbwithparts, gbc, gpc, est, gss, seqid, acc, ft. Formerly, but not currently, ‘genbank’ was a synonym for ‘gb’ and ‘genpept’ was a synonym for ‘gp’: however, these synonyms no longer work and need to be fixed if you encounter them in old code. For more information check NCBI’s format documentation.

Note that there are two separate concepts: rettype and retmode. rettype controls what kind of data you’ll get, and retmode controls how the data will be formatted.

For example:

    >>> from cogent.db.ncbi import EFetch
    >>> ef = EFetch(id='459567', rettype='fasta', retmode='text')
    >>> lines = ef.read().splitlines()
    >>> for line in lines:
    ...     print line[:40]
    ...
    >gi|459567|dbj|D28543.1|HPCNS5PC Hepatit
    GAGCACGACATCTACCAATGTTGCCAACTGAACCCAGAGG
    GGCTTTACCTTGGTGGTCCCATGTTTAACTCGCGAGGTCA
    CGGGGTTCTTCCAACCAGCATGGGCAATACCCTCACATGT
    GCAGGCCTCACCAATTCTGACATGTTGGTTTGCGGAGATG
    TC

    >>> ef = EFetch(id='459567', rettype='fasta', retmode='html')
    >>> lines = ef.read().splitlines()
    >>> for line in lines:
    ...     print line[:40]
    ...
    >gi|459567|dbj|D28543.1|HPCNS5PC Hepatit
    GAGCACGACATCTACCAATGTTGCCAACTGAACCCAGAGG
    GGCTTTACCTTGGTGGTCCCATGTTTAACTCGCGAGGTCA
    CGGGGTTCTTCCAACCAGCATGGGCAATACCCTCACATGT
    GCAGGCCTCACCAATTCTGACATGTTGGTTTGCGGAGATG
    TC

    >>> ef = EFetch(id='459567', rettype='fasta', retmode='xml')
    >>> lines = ef.read().splitlines()
    >>> for line in lines:
    ...     print line[:40]
    ...
    <?xml version="1.0"?>
     <!DOCTYPE TSeqSet PUBLIC "-//NCBI//NCBI
     <TSeqSet>
    <TSeq>
      <TSeq_seqtype value="nucleotide"/>
      <TSeq_gi>459567</TSeq_gi>
      <TSeq_accver>D28543.1</TSeq_accver>
      <TSeq_taxid>11103</TSeq_taxid>
      <TSeq_orgname>Hepatitis C virus</TSeq_
      <TSeq_defline>Hepatitis C virus gene f
      <TSeq_length>282</TSeq_length>
      <TSeq_sequence>GAGCACGACATCTACCAATGTTG
    </TSeq>

    </TSeqSet>

You’ll notice that the second case is some funny-looking html. Thanks, NCBI! This is not our fault, please don’t file a bug report. To figure out whether something is actually surprising behavior at NCBI, you can always capture the command-line and run it in a web browser. You can do this by calling str() on the `ef`, or by printing it. For example:

    >>> print ef
    http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?retmax=100&retmod...

If you paste the resulting string into your web browser and you get the same incorrect result that you get using PyCogent, you know that you should direct your support requests NCBI’s way. If you want to use your own email address instead of leaving it as the default (the module developer), you can do that just by passing it in as a parameter. For example, in the unlikely event that I want NCBI to contact me instead of Mike if something goes wrong with my script, I can achieve that as follows:

    >>> ef = EFetch(id='459567', rettype='fasta', retmode='xml', email='rob@spot.colorado.edu')
    >>> print ef
    http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?retmax=100&retmod...

You can also select multiple ids (pass in as comma-delimited list):

    >>> ef = EFetch(id='459567,459568', rettype='brief')
    >>> ef.read()
    'D28543 Hepatitis C virus... [gi:459567]\nBAA05896 NS5 protein [Hepa... [gi:459568]'

#### Retrieving GenPept files from NCBI via Eutils

We query for just one accession to illustrate the process. A more general query can be executed by replacing `'BAB52044` with `'"lysyl tRNA-synthetase"[ti] AND bacteria[orgn]'` in the snippet below.

    >>> from cogent.db.ncbi import EUtils
    >>> e = EUtils(numseqs=100, db='protein', rettype='gp')
    >>> result = e['BAB52044']
    >>> print result.read()
    LOCUS       BAB52044                 548 aa            linear   BCT 16-MAY-2009
    DEFINITION  lysyl tRNA synthetase [Mesorhizobium loti MAFF303099].
    ACCESSION   BAB52044
    VERSION     BAB52044.1  GI:14025444
    DBSOURCE    accession BA000012.4
    KEYWORDS    .
    SOURCE      Mesorhizobium loti MAFF303099...

#### Retrieving and parsing GenBank entries

    >>> from cogent.parse.genbank import RichGenbankParser
    >>> from cogent.db.ncbi import EUtils
    >>> e = EUtils(numseqs=100, db='protein', rettype='gp')
    >>> result = e['"lysyl tRNA-synthetase"[ti] AND bacteria[orgn]']
    >>> parser = RichGenbankParser(result.readlines())
    >>> gb = [(accession, seq) for accession, seq in parser]

Printing the resulting list (`gb`) will generate output like:

    [('AAA83071', Sequence(MSEQHAQ... 505)), ('ACS40931', ...

#### Parsing in more detail: a single GenBank entry

    >>> from cogent.db.ncbi import EUtils
    >>> from cogent.parse.genbank import RichGenbankParser
    >>> e = EUtils(db="nucleotide", rettype="gb")
    >>> record = e['154102'].readlines()
    >>> parser = RichGenbankParser(record)
    >>> accession, seq = [record for record in parser][0]
    >>> accession
    'STYHEMAPRF'
    >>> type(seq)
    <class 'cogent.core.sequence.DnaSequence'>
    >>> def gene_and_cds(f):
    ...     return f['type'] == 'CDS' and 'gene' in f
    ...
    >>> cds_features = [f for f in seq.Info.features if gene_and_cds(f)]
    >>> for cds_feature in cds_features:
    ...     print cds_feature['gene'], cds_feature['location']
    ...
    ['hemA'] 732..1988
    ['prfA'] 2029..3111

#### Retrieving a bacterial genome file

To obtain a full bacterial genome, run the following to get the complete Salmonella typhimurium genome sequence (Genbank) file. (For this documentation, we include a partial file for illustration purposes.)

    from cogent.db.ncbi import EUtils
    e = EUtils(db="nucleotide", rettype="gb")
    outfile = open('data/ST.genome.gb','w')
    outfile.write(e['AE006468'].read())
    outfile.close()

For larger files, you might want to dump them directly into a file on your hard drive rather than reading them into memory first, e.g.

    e.filename='ST.genome.gb'
    f = e['AE006468']

dumps the result into the file directly, and returns you a handle to the open file where you can read the result, get the path, or do any of the other standard file operations. Now do the analysis:

    >>> from cogent.parse.genbank import RichGenbankParser
    >>> infile = open('data/ST_genome_part.gb', 'r')
    >>> parser = RichGenbankParser(infile)
    >>> accession, seq = [record for record in parser][0]
    >>> gene_and_cds = lambda f: f['type'] == 'CDS' and 'gene' in f
    >>> gene_name = lambda f: f['gene']
    >>> all_cds = [f for f in seq.Info.features if gene_and_cds(f)]
    >>> for cds in sorted(all_cds, key=gene_name):
    ...     print cds['gene'][0].ljust(6),
    ...     print cds['protein_id'], cds['location']
    ...
    mog    ['AAL18972.1'] 8729..9319
    talB   ['AAL18971.1'] 7665..8618
    thrA   ['AAL18966.1'] 337..2799
    thrB   ['AAL18967.1'] 2801..3730
    thrC   ['AAL18968.1'] 3734..5020
    thrL   ['AAL18965.1'] 190..255
    yaaA   ['AAL18969.1'] complement(5114..5887)
    yaaH   ['AAL18973.1'] complement(9376..9942)
    yaaJ   ['AAL18970.1'] complement(5966..7396)

The EUtils modules are generic, so additional databases like OMIM can be accessed using similar mechanisms.

#### Retrieving PubMed abstracts from NCBI via EUtils

    >>> from cogent.db.ncbi import EUtils
    >>> e = EUtils(db='pubmed',rettype='brief')
    >>> result = e['Simon Easteal AND Von Bing Yap'].read()
    >>> print result

    1: Yap VB et al. Estimates of the effect of na...[PMID: 19815689]

    2: Schranz HW et al. Pathological rate matrices: f...[PMID: 19099591]

#### Retrieving PubMed abstracts via PMID

    >>> from cogent.db.ncbi import EUtils
    >>> e = EUtils(db='pubmed',rettype='abstract')
    >>> result = e['14983078'].read()

## KEGG

### Complete genomes

To be written.

### Orthologs

To be written.

### Functional assignments

To be written.

### Pathway assignments

To be written.

## Ensembl

Note that much more extensive documentation is available in Querying Ensembl.

### Connecting

Ensembl provides access to their MySQL databases directly or users can download and run those databases on a local machine. To use the Ensembl’s UK servers for running queries, nothing special needs to be done as this is the default setting for PyCogent’s `ensembl` module. To use a different Ensembl installation, you create an account instance:

    >>> from cogent.db.ensembl import HostAccount
    >>> account = HostAccount('fastcomputer.topuni.edu', 'username',
    ...                       'canthackthis')

To specify a specific port to connect to MySQL on:

    >>> from cogent.db.ensembl import HostAccount
    >>> account = HostAccount('anensembl.server.edu', 'someuser',
    ...                       'somepass', port=3306)

### Species to be queried

To see what existing species are available

    >>> from cogent.db.ensembl import Species
    >>> print Species
    ================================================================================
           Common Name                   Species Name              Ensembl Db Prefix
    --------------------------------------------------------------------------------
             A.aegypti                  Aedes aegypti                  aedes_aegypti
            A.clavatus           Aspergillus clavatus           aspergillus_clavatus...

If Ensembl has added a new species which is not yet included in `Species`, you can add it yourself.

    >>> Species.amendSpecies('A latinname', 'a common name')

You can get the common name for a species

    >>> Species.getCommonName('Procavia capensis')
    'Rock hyrax'

and the Ensembl database name prefix which will be used for all databases for this species.

    >>> Species.getEnsemblDbPrefix('Procavia capensis')
    'procavia_capensis'

Species common names are used to construct attributes on PyCogent `Compara` instances). You can get the name that will be using the `getComparaName` method. For species with a real common name

    >>> Species.getComparaName('Procavia capensis')
    'RockHyrax'

or with a shortened species name

    >>> Species.getComparaName('Caenorhabditis remanei')
    'Cremanei'

### Get genomic features

#### Find a gene by gene symbol

We query for the BRCA2 gene for humans.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> print human
    Genome(Species='Homo sapiens'; Release='67')
    >>> genes = human.getGenesMatching(Symbol='BRCA2')
    >>> for gene in genes:
    ...     if gene.Symbol == 'BRCA2':
    ...         print gene
    ...         break
    Gene(Species='Homo sapiens'; BioType='protein_coding'; Description='breast cancer 2,...'; StableId='ENSG00000139618'; Status='KNOWN'; Symbol='BRCA2')

#### Find a gene by Ensembl Stable ID

We use the stable ID for BRCA2.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> gene = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> print gene
    Gene(Species='Homo sapiens'; BioType='protein_coding'; Description='breast cancer 2,...'; StableId='ENSG00000139618'; Status='KNOWN'; Symbol='BRCA2')

#### Find genes matching a description

We look for breast cancer related genes that are estrogen induced.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> genes = human.getGenesMatching(Description='breast cancer anti-estrogen')
    >>> for gene in genes:
    ...     print gene
    Gene(Species='Homo sapiens'; BioType='lincRNA'; Description='breast cancer anti-estrogen...'; StableId='ENSG00000262117'; Status='NOVEL'; Symbol='BCAR4')...

We can also require that an exact (case insensitive) match to the word(s) occurs within the description by setting `like=False`.

    >>> genes = human.getGenesMatching(Description='breast cancer anti-estrogen',
    ...                                  like=False)
    >>> for gene in genes:
    ...     print gene
    Gene(Species='Homo sapiens'; BioType='lincRNA'; Description='breast cancer anti-estrogen...'; StableId='ENSG00000262117'; Status='NOVEL'; Symbol='BCAR4')...

#### Get canonical transcript for a gene

We get the canonical transcripts for BRCA2.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> transcript = brca2.CanonicalTranscript
    >>> print transcript
    Transcript(Species='Homo sapiens'; CoordName='13'; Start=32889610; End=32973347; length=83737; Strand='+')

#### Get the CDS for a transcript

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> transcript = brca2.CanonicalTranscript
    >>> cds = transcript.Cds
    >>> print type(cds)
    <class 'cogent.core.sequence.DnaSequence'>
    >>> print cds
    ATGCCTATTGGATCCAAAGAGAGGCCA...

#### Look at all transcripts for a gene

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> for transcript in brca2.Transcripts:
    ...     print transcript
    Transcript(Species='Homo sapiens'; CoordName='13'; Start=32889610; End=32973347; length=83737; Strand='+')
    Transcript(Species='Homo sapiens'; CoordName='13'; Start=32889641; End=32907428; length=17787; Strand='+')...

#### Get the first exon for a transcript

We show just for the canonical transcript.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> print brca2.CanonicalTranscript.Exons[0]
    Exon(StableId=ENSE00001184784, Rank=1)

#### Get the introns for a transcript

We show just for the canonical transcript.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> for intron in brca2.CanonicalTranscript.Introns:
    ...     print intron
    Intron(TranscriptId=ENST00000380152, Rank=1)
    Intron(TranscriptId=ENST00000380152, Rank=2)
    Intron(TranscriptId=ENST00000380152, Rank=3)...

#### Inspect the genomic coordinate for a feature

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> print brca2.Location.CoordName
    13
    >>> print brca2.Location.Start
    32889610
    >>> print brca2.Location.Strand
    1

#### Get repeat elements in a genomic interval

We query the genome for repeats within a specific coordinate range on chromosome 13.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> repeats = human.getFeatures(CoordName='13', Start=32879610, End=32889610, feature_types='repeat')
    >>> for repeat in repeats:
    ...     print repeat.RepeatClass
    ...     print repeat
    ...     break
    SINE/Alu
    Repeat(CoordName='13'; Start=32879362; End=32879662; length=300; Strand='-', Score=2479.0)

#### Get CpG island elements in a genomic interval

We query the genome for CpG islands within a specific coordinate range on chromosome 11.

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> islands = human.getFeatures(CoordName='11', Start=2150341, End=2170833, feature_types='cpg')
    >>> for island in islands:
    ...     print island
    ...     break
    CpGisland(CoordName='11'; Start=2158951; End=2162484; length=3533; Strand='-', Score=3254.0)

### Get SNPs

#### For a gene

We find the genetic variants for the canonical transcript of BRCA2.

Note

The output is significantly truncated!

    >>> from cogent.db.ensembl import Genome
    >>> human = Genome('human', Release=67, account=account)
    >>> brca2 = human.getGeneByStableId(StableId='ENSG00000139618')
    >>> transcript = brca2.CanonicalTranscript
    >>> print transcript.Variants
    (<cogent.db.ensembl.region.Variation object at ...
    >>> for variant in transcript.Variants:
    ...     print variant
    ...     break
    Variation(Symbol='rs55880202'; Effect=['2KB_upstream_variant', '5_prime_UTR_variant', '5KB_upstream_variant']; Alleles='C/T')...

#### Get a single SNP

We get a single SNP and print it’s allele frequencies.

    >>> snp = list(human.getVariation(Symbol='rs34213141'))[0]
    >>> print snp.AlleleFreqs
    =============================
    allele      freq    sample_id
    -----------------------------
         A    0.0303          933
         G    0.9697          933
         G    1.0000        11208
         G    1.0000        11519
         A    0.0110       113559
         G    0.9889       113559...

### What alignment types available

We create a `Compara` instance for human, chimpanzee and macaque.

    >>> from cogent.db.ensembl import Compara
    >>> compara = Compara(['human', 'chimp', 'macaque'], Release=67,
    ...                  account=account)
    >>> print compara.method_species_links
    Align Methods/Clades
    ===================================================================================================================
    method_link_species_set_id  method_link_id  species_set_id      align_method                            align_clade
    -------------------------------------------------------------------------------------------------------------------
                           580              10           34468             PECAN           19 amniota vertebrates Pecan
                           548              13           34115               EPO                         6 primates EPO
                           578              13           34466               EPO               12 eutherian mammals EPO
                           582              14           34697  EPO_LOW_COVERAGE  35 eutherian mammals EPO_LOW_COVERAGE
    -------------------------------------------------------------------------------------------------------------------

### Get genomic alignment for a gene region

We first get the syntenic region corresponding to human gene BRCA2.

    >>> from cogent.db.ensembl import Compara
    >>> compara = Compara(['human', 'chimp', 'macaque'], Release=67,
    ...                  account=account)
    >>> human_brca2 = compara.Human.getGeneByStableId(StableId='ENSG00000139618')
    >>> regions = compara.getSyntenicRegions(region=human_brca2, align_method='EPO', align_clade='primates')
    >>> for region in regions:
    ...     print region
    SyntenicRegions:
      Coordinate(Human,chro...,13,32889610-32973805,1)
      Coordinate(Macaque,chro...,17,11686607-11778803,1)
      Coordinate(Chimp,chro...,13,31957346-32040817,1)...

We then get a cogent `Alignment` object, requesting that sequences be annotated for gene spans.

    >>> aln = region.getAlignment(feature_types='gene')
    >>> print repr(aln)
    3 x 98805 dna alignment: Homo sapiens:chromosome:13:3288...

### Parsing syntenic regions

Not all regions in a given genome have a syntenic alignment, and some have more than one alignment. To illustrate these cases, we can consider an alignment between mouse and human, using the `PECAN` alignment method in the vertebrates clade:

    >>> species = ["mouse", "human"]
    >>> compara = Compara(species, Release=67, account=account)
    >>> clade = "vertebrates"
    >>> chrom, start, end, strand = "X", 165754928, 165755079, "-"
    >>> regions = compara.getSyntenicRegions(Species="mouse", CoordName=chrom,
    ...                                      Start=start, End=end, align_method="PECAN",
    ...                                      align_clade=clade, Strand=strand)
    >>> aligned_pairs = [r for r in regions]
    >>> alignment = aligned_pairs[0]
    >>> aligned_regions = [m for m in alignment.Members
    ...                    if m.Region is not None]
    >>> source_region, target_region = aligned_regions
    >>> print source_region.Location.CoordName, source_region.Location.Start, source_region.Location.End
    X 165754928 165755079
    >>> print target_region.Location.CoordName, target_region.Location.Start, target_region.Location.End
    X 11132954 11133105

Note

We took the aligned regions from the `regions` generator and put them in a list for convenience.

If there are no regions returned (i.e. `num_pairs` is zero), then no alignment could be found. In the case of the above region, an exon in the Hccs gene, there is only one alignment. We then accessed the coordinates of the alignment using the `Members` attribute of the region. Each element of `aligned_regions` is a `SyntenicRegion` instance, whose coordinates can be pulled from the `Location` attribute.

This example shows that mouse region `X:165754928-165755079` aligns only to human region `X:11132954-11133105`.

Note

Sometimes, the genomic coordinates given to `getSyntenicRegions` will contain multiple alignments between the pair of genomes, in which case two or more regions will be returned in `aligned_pairs`.

### Getting related genes

#### What gene relationships are available

    >>> from cogent.db.ensembl import Compara
    >>> compara = Compara(['human', 'chimp', 'macaque'], Release=67,
    ...                  account=account)
    >>> print compara.getDistinct('relationship')
    [u'ortholog_one2many', u'contiguous_gene_split', u'ortholog_one2one',...

#### Get one-to-one orthologs

We get the one-to-one orthologs for BRCA2.

    >>> from cogent.db.ensembl import Compara
    >>> compara = Compara(['human', 'chimp', 'macaque'], Release=67,
    ...                  account=account)
    >>> orthologs = compara.getRelatedGenes(StableId='ENSG00000139618',
    ...                  Relationship='ortholog_one2one')
    >>> print orthologs
    RelatedGenes:
     Relationships=ortholog_one2one
      Gene(Species='Macaca mulatta'; BioType='protein_coding'; Description=...

We iterate over the related members.

    >>> for ortholog in orthologs.Members:
    ...     print ortholog
    Gene(Species='Macaca mulatta'; BioType='protein_coding'; Description=...

We get statistics on the ortholog CDS lengths.

    >>> print orthologs.getMaxCdsLengths()
    [10008, 10257, 10257]

We get the sequences as a sequence collection, with annotations for gene.

    >>> seqs = orthologs.getSeqCollection(feature_types='gene')

#### Get CDS for all one-to-one orthologs

We sample all one-to-one orthologs for a group of species, generating a FASTA formatted string that can be written to file. We check all species have an ortholog and that all are translatable.

    >>> from cogent.core.alphabet import AlphabetError
    >>> common_names = ["mouse", "rat", "human", "opossum"]
    >>> latin_names = set([Species.getSpeciesName(n) for n in common_names])
    >>> latin_to_common = dict(zip(latin_names, common_names))
    >>> compara = Compara(common_names, Release=67, account=account)
    >>> for gene in compara.Human.getGenesMatching(BioType='protein_coding'):
    ...     orthologs = compara.getRelatedGenes(gene,
    ...                                  Relationship='ortholog_one2one')
    ...     # make sure all species represented
    ...     if orthologs is None or orthologs.getSpeciesSet() != latin_names:
    ...         continue
    ...     seqs = []
    ...     for m in orthologs.Members:
    ...         try: # if sequence can't be translated, we ignore it
    ...             # get the CDS without the ending stop
    ...             seq = m.CanonicalTranscript.Cds.withoutTerminalStopCodon()
    ...             # make the sequence name
    ...             seq.Name = '%s:%s:%s' % \
    ...         (latin_to_common[m.genome.Species], m.StableId, m.Location)
    ...             aa = seq.getTranslation()
    ...             seqs += [seq]
    ...         except (AlphabetError, AssertionError):
    ...             seqs = [] # exclude this gene
    ...             break
    ...     if len(seqs) == len(common_names):
    ...         fasta = '\n'.join(s.toFasta() for s in seqs)
    ...         break

#### Get within species paralogs

    >>> paralogs = compara.getRelatedGenes(StableId='ENSG00000164032',
    ...             Relationship='within_species_paralog')
    >>> print paralogs
    RelatedGenes:
     Relationships=within_species_paralog
      Gene(Species='Homo sapiens'; BioType='protein_coding'; Description='H2A histone...

## PDB

### For structures

The PDB module is very simple and basically gets a pdb coordinates file by accession number. Searches, etc. are not currently implemented because it’s easier to get the pdb ids from NCBI than to scrape PDB’s html results format.

    >>> from cogent.db.pdb import Pdb
    >>> p = Pdb()
    >>> result = p['3L0U']

returns a handle to a file containing the PDB coordinates (that you can, for example, pass to the PDB parser in a fashion analogous to how you pass the GenBank record above to the RichGenbankParser). See the pdb parser documentation for more info. To send results directly to a file, you can use the retrieve() method of the Pdb object.

## Rfam

### For RNA secondary structures, alignments, functions

To be written.

## GoldenPath (not yet implemented)

To be written.

### Whole-genome alignments, orthologs, annotation tracks

To be written.