run_geo_reuse_collection.py

import os
import time
import re
import nose
from nose.tools import assert_equals
from collections import defaultdict
import EUtils
from EUtils import HistoryClient, ThinClient
import dataset
import datasources
from datasources import pubmedcentral
from datasources import pubmed
from datasources import affiliation
from datasources import geo
from datasources import urlopener
import utils
from utils.cache import TimedCache
import pickle

EMAIL_CONTACT = "hpiwowar@gmail.com"
VERBOSE = False

#base_query = """(GEO[text] OR omnibus[text]) NOT "pmc gds"[filter]"""
#base_query = """(GEO[text] OR omnibus[text]) NOT "pmc gds"[filter] AND ("1900"[PubDate] : "2009"[PubDate])"""
base_query_reuse = """("1900"[PubDate] : "2009"[PubDate]) NOT "pmc gds"[filter]"""
base_query_submit = """("1900"[PubDate] : "2009"[PubDate]) AND "pmc gds"[filter]"""

#url_for_gse = """http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gds&term=GSE[ETYP]&retmax=10000&usehistory=n"""
#url_for_gds = """http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gds&term=GDS[ETYP]&retmax=10000&usehistory=n"""


def test_get_accession_in_pmc_fulltext():
    response = get_accession_in_pmc_fulltext("GSE", "200008514", base_query_reuse)
    assert_equals(response, ['2785812', '2620272'])

    response = get_accession_in_pmc_fulltext("GSE", "200008478", base_query_reuse)
    assert_equals(response, ['2223695'])

    response = get_accession_in_pmc_fulltext("GSE", "200007778", base_query_reuse)
    assert_equals(response, [])

    response = get_accession_in_pmc_fulltext("GDS", "2643", base_query_reuse)
    assert_equals(response, ['2715883', '2781753'])
    
def test_get_accession_variants():
    response = get_accession_variants("GSE", "200008514")
    assert_equals(response, ['GSE200008514', '"GSE 200008514"', 'GSE8514', '"GSE 8514"'])
    
def get_accession_variants(id_type, id):
    variants = []
    variants.append(id_type + id)
    variants.append('"' + id_type + ' ' + id  + '"')
    id_stripped = geo.get_stripped_accession(id)
    if (id_stripped != id):
        variants.append(id_type + id_stripped)
        variants.append('"' + id_type + ' ' + id_stripped + '"')
    return(variants)
        
def get_accession_in_pmc_fulltext(id_type, id, pmc_query):
    pmc_ids = []
    accession_variants = get_accession_variants(id_type, id)
    for variant in accession_variants:
        query = variant + "[text] AND " + pmc_query
        pmc_ids += pubmedcentral.search(query)        
    return(pmc_ids)
    
def test_get_authors_and_submittors_from_accession():
    response = get_authors_and_submittors_from_accession([u'17911395'], [u'Ye P', u'Rainey W'])
    assert_equals(response, ['Mariniello', 'Shibata', 'Ye', 'Rainey', 'Mantero'])

def get_author_last_names(pmids):
    last_names_list = []
    for pmid in pmids:
        authors = pubmed.authors([pmid])
        authors_list = authors[0].split(";")
        last_names = [author.split(" ")[0] for author in authors_list]
        last_names_list += last_names
    return(last_names_list)
    
def get_authors_and_submittors_from_accession(pmids, contributors=None):
    last_names_list = get_author_last_names(pmids)
    if contributors:
        contributor_last_names = [author.split(" ")[0] for author in contributors]
        last_names_list += contributor_last_names
    last_names_set = set(last_names_list)
    return(list(last_names_set))
    
def test_get_dict_submit_reuse_of_accession_in_pmc_fulltext():
    response = get_dict_submit_reuse_of_accession_in_pmc_fulltext("GSE", ["200008514", "200008478"], base_query_reuse)
    assert_equals(response.items()[0:2], [(('GSE200008514', '2785812'), ('200008514', [u'17911395'], ['2785812', '2620272'], [u'19917117'], ['Mariniello', 'Shibata', 'Ye', 'Rainey', 'Mantero'], ['Kroll', 'Barkema', 'Carlon'], [], ['Algorithms Mol Biol'], ['2009'], [], [], [])), (('GSE200008478', '2223695'), ('200008478', [u'17977147', u'17951393'], ['2223695'], [u'17993534'], ['Moser', 'Fischer', 'Friberg', 'Pessi', 'Lindemann', 'Hauser', 'Rehrauer', 'Hennecke', 'Ahrens'], ['Hacker', 'Aktas', 'Narberhaus', 'Sohlenkamp', 'Geiger'], [], ['J Bacteriol'], ['2008'], [], [], ['deposited', 'accessed', 'published']))] )
    
def test_get_author_intersect_submit_reuse():
    response = get_author_intersect_submit_reuse(['Pessi', 'Ahrens', 'Rehrauer', 'Lindemann', 'Hauser', 'Fischer', 'Hennecke', 'Lindemann', 'Moser', 'Pessi', 'Hauser', 'Friberg', 'Hennecke', 'Fischer', u'Pessi', u'Ahrens', u'Rehrauer', u'Lindemann', u'Hauser', u'Fischer', u'Hennecke'], ['Hacker', 'Sohlenkamp', 'Pessi', 'Aktas', 'Geiger', 'Narberhaus'])
    assert_equals(response, ['Pessi'])
    
def get_author_intersect_submit_reuse(submit_authors, reuse_authors):
    intersect = multi_intersection([submit_authors, reuse_authors])
    return(intersect)

meshes = """"Algorithms"[mesh]
"Databases, Genetic"[mesh]
"Gene Expression Profiling/methods"[mesh]
"Computational Biology/methods"[mesh]
"Oligonucleotide Array Sequence Analysis/methods"[mesh]
"Genomics/methods"[mesh]
"Reproducibility of Results"[mesh]
"Software"[mesh]
"Computer Simulation"[mesh]
"Internet"[mesh]
"Data Interpretation, Statistical"[mesh]""".split("\n")

metaquery = "(meta-analysis [pt] OR meta-analysis [tw] OR metanalysis [tw]) OR meta-analysis [mh])"

words = """submitted
deposited
user*
public
accessed
downloaded
published""".split("\n")
        
def get_dict_submit_reuse_of_accession_in_pmc_fulltext(id_type, ids, pmc_query):
    response_dict = defaultdict(str)
    num_ids = len(ids)
    id_counter = 0
    geo_instance = geo.GEO()
    for accession in ids:
        id_counter += 1
        pmc_reuse_pmcids = get_accession_in_pmc_fulltext(id_type, accession, pmc_query)        
        if not pmc_reuse_pmcids:
            continue
        stripped_accession = geo.get_stripped_accession(accession)     
        try:       
            submit_pmids = geo_instance.pmids(id_type + stripped_accession)
        except Exception:
            continue
        submit_contributors = geo_instance.contributors(id_type + stripped_accession)
        submit_authors = get_authors_and_submittors_from_accession(submit_pmids, submit_contributors)

        for reuse_pmcid in pmc_reuse_pmcids:
            reuse_pmids      = pubmedcentral.pmcids_to_pmids(reuse_pmcid)
            reuse_authors   = get_authors_and_submittors_from_accession(reuse_pmids)
            submit_affiliation = affiliation.institution(submit_pmids) if submit_pmids else []
            reuse_affiliation = affiliation.institution(reuse_pmids)
            intersect       = get_author_intersect_submit_reuse(submit_authors, reuse_authors)
            journal         = pubmed.journal(reuse_pmids)
            year            = pubmed.year_published(reuse_pmids)
            medline_status  = pubmed.medline_status(reuse_pmids)
            metaanal        = pubmed.filter_pmids(reuse_pmids, metaquery)
            biolink_filter    = pubmedcentral.filter_pmcids([reuse_pmcid], '(geo OR omnibus)  AND microarray  AND "gene expression" AND accession NOT (databases OR user OR users  OR (public AND accessed) OR (downloaded AND published))')
            basic_filter  = pubmedcentral.filter_pmcids([reuse_pmcid], '"gene expression omnibus" AND (submitted OR deposited)')
            mesh_filters    = [term for term in meshes if pubmed.filter_pmids(reuse_pmids, term)]
            word_filters    = [term for term in words if pubmedcentral.filter_pmcids([reuse_pmcid], term + "[text]")]
        
            response_dict[(id_type+stripped_accession, reuse_pmcid)] = (id_type+stripped_accession, submit_pmids, reuse_pmcid, reuse_pmids, submit_authors, reuse_authors, submit_affiliation, reuse_affiliation, intersect, journal, year, medline_status, metaanal, biolink_filter, basic_filter, mesh_filters, word_filters)
            print response_dict[(id_type+stripped_accession, reuse_pmcid)]
            print id_counter, "of", num_ids, ":", stripped_accession, "--", (submit_pmids), "; ", len(pmc_reuse_pmcids)
    return(response_dict)
    
def authors_in_common_from_pmids():
    pmids = ["20349403", "18998887", "18767901"]
    response = authors_in_common_from_pmids(pmids)
    assert_equals(response, ["Piwowar"])

def multi_intersection(xs):
    inter = reduce(set.intersection, [set(x) for x in xs])
    return list(inter)
    
def authors_in_common_from_pmids(pmids):
    last_names_list = []
    for pmid in pmids:
        authors = pubmed.authors([pmid])
        authors_list = authors[0].split(";")
        last_names = [author.split(" ")[0] for author in authors_list]
        last_names_list.append(last_names)
    authors_intersection = multi_intersection(last_names_list)
    return(authors_intersection)
    
def get_from_query_gds_in_pmc_fulltext_dict(id_type, pmc_query, geo_year=None):  
    ids = geo.get_ids_by_year(id_type, geo_year)
    (response_dict) = get_dict_submit_reuse_of_accession_in_pmc_fulltext(id_type, ids, pmc_query)  
    return(response_dict)


def print_accession_pmcids(prefix, accession_dict):
    for accession in accession_dict:
        pmcid_list = accession_dict[accession]
        if pmcid_list:
            for pmcid in pmcid_list:
                print "%s\t%s%s\t%s" %(prefix, prefix, accession, pmcid)
        else:
            print "%s\t%s%s\t%s" %(prefix, prefix, accession, "")

#print_accession_pmcids("GDS", gds_dict)
#print_accession_pmcids("GSE", gse_dict)

def test_estimate_pmc_coverage():
    response = estimate_pmc_coverage("arrayexpress[title]")
    assert_equals(response, (6, 14, 0.42857142857142855))

    response = estimate_pmc_coverage('"gene expression profiling"[mesh]')
    assert_equals(response, (11038, 47384, 0.23294783049130507))

def test_estimate_pmc_coverage_given_years():
    response = estimate_pmc_coverage('"gene expression profiling"[mesh]', "2007", "2009")
    assert_equals(response, (6311, 21569, 0.29259585516250175))
    
def estimate_pmc_coverage(query, start_year="1800", end_year="3000"):
    pubmed_query = query + ' AND ("' + start_year + '"[pdat] : "' + end_year + '"[pdat])'
    pubmed_ids = pubmed.search(pubmed_query)
    num_pubmed = len(pubmed_ids)

    pmc_query = query + ' AND ("' + start_year + '"[PubDate] : "' + end_year + '"[PubDate])'
    pmc_ids = pubmedcentral.search(pmc_query)
    num_pmc = len(pmc_ids)
    
    ratio = num_pmc / (num_pubmed + 0.0)
    
    return(num_pmc, num_pubmed, ratio)

def run_stats():
    if False:
        (num_pmc, num_pubmed, ratio) = estimate_pmc_coverage('"gene expression profiling"[mesh]', "2007", "2009")

    geo_year = "2007"
    id_types = ["GDS", "GSE"]
    response_dict = {}
    for id_type in id_types:
        response_dict[(id_type, geo_year)] = get_from_query_gds_in_pmc_fulltext_dict(id_type, base_query_reuse, geo_year)
        pkl_file = open("scienceplot/results/" + id_type + "_dict" + geo_year + ".pkl", "wb")
        pickle.dump(response_dict[(id_type, geo_year)], pkl_file)
        pkl_file.close()
    
        fh = open("scienceplot/results/" + id_type + geo_year + ".csv", "w")
        header = "id_type+stripped_accession, submit_pmids, reuse_pmcid, reuse_pmids, submit_authors, reuse_authors, submit_affiliation, reuse_affiliation, intersect, journal, year, medline_status, metaanal, biolink_filter, basic_filter, mesh_filters, word_filters"
        fh.write(header + "\n")
        dataset.csv_write_to_file(fh, response_dict[(id_type, geo_year)].values())
        fh.close()

    return(response_dict)