功能富集作为生信分析的常规方法,做过分析的人都很熟悉。功能富集的软件也是相当的丰富。其中,DAVID作为功能富集的老牌网站认可度还是很高的,对于不熟悉代码的人来说,这样可视化的网站简直是不可多得的选择,点点鼠标即可获得想要的富集结果,可以说是不费吹灰之力。

虽说网页可以满足富集分析的整体任务,但个人还是喜欢用代码的方式来解决,因为使用代码的方式有几个好处:可以将分析过程保存下来方便复盘回溯;可重现性及反复利用;可以在此基础上很方便地继续进行个性化修改。
Some Limitations:
- A job with more than 3000 genes to generate gene or term cluster report will not be handled by DAVID due to resource limit.
- No more than 200 jobs in a day from one user or computer.
- DAVID Team reserves right to suspend any improper uses of the web service without notice.
想要使用DAVID Web Service功能,需要先到下面的网址完成注册:
https://davidbioinformatics.nih.gov/webservice/register.htm
from suds.client import Client
import pandas as pd
df = pd.read_csv('demolist.txt', usecols=[0], delimiter='\t', index_col=False, names=['AFFYMETRIX_3PRIME_IVT_ID'])
inputIds = ",".join(df['AFFYMETRIX_3PRIME_IVT_ID'].astype(str).unique().tolist())
client = Client("https://davidbioinformatics.nih.gov/webservice/services/DAVIDWebService?wsdl")
client.wsdl.services[0].setlocation('https://davidbioinformatics.nih.gov/webservice/services/DAVIDWebService.DAVIDWebServiceHttpSoap11Endpoint')
client.service.authenticate('YourRegisteredEmail@your.org')
client.service.addList(inputIds, 'AFFYMETRIX_3PRIME_IVT_ID', 'genelist', 0)
chart = client.service.getChartReport(0.1, 2)
df = pd.DataFrame([{"Category": x.categoryName, "Term": x.termName, "Count": x.listHits, "Percent": x.percent, "PValue": x.ease, "Genes": x.geneIds, "Totals":x.listTotals, "popHits": x.popHits, "foldEnrichment": x.foldEnrichment, "Benjamini": x.benjamini, "FDR": x.afdr} for x in chart])
df.head()
Category Term Count Percent PValue Genes Totals popHits foldEnrichment Benjamini FDR
0 GOTERM_BP_DIRECT GO:0061844~antimicrobial humoral immune respon... 10 6.944444 4.789319e-07 35367_AT, 38604_AT, 37454_AT, 1403_S_AT, 34623... 134 137 10.443403 0.000611 0.060537
1 GOTERM_CC_DIRECT GO:0005576~extracellular region 45 31.250000 9.956898e-07 1645_AT, 34623_AT, 1372_AT, 35091_AT, 36766_AT... 137 3186 2.092767 0.000250 0.023996
2 UP_KW_CELLULAR_COMPONENT KW-0964~Secreted 36 25.000000 1.369197e-06 1645_AT, 34623_AT, 1372_AT, 35091_AT, 37898_R_... 125 2237 2.346357 0.000052 0.005203
3 GOTERM_BP_DIRECT GO:0050829~defense response to Gram-negative b... 8 5.555556 8.536393e-06 38604_AT, 34623_AT, 679_AT, 33963_AT, 33979_AT... 134 105 10.900924 0.005446 0.539500
4 GOTERM_BP_DIRECT GO:0019731~antibacterial humoral response 7 4.861111 1.963736e-05 38604_AT, 34623_AT, 679_AT, 33963_AT, 33979_AT... 134 80 12.519030 0.008352 0.827387
上面的代码就是做富集分析的过程,把这些代码复制一下然后往python解释器一扔结果就出来,这可比网页操作快捷多了。如果想基于结果继续搞点啥,直接搞就行。
SSLError: [SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1000)
如果创建网络连接时出现类型上面的错误,可用下面的临时解决方法:
import ssl
ssl._create_default_https_context = ssl._create_unverified_context
如果想自定义背景基因集可用下面的方式:
client.service.addList(bgenes, 'AFFYMETRIX_3PRIME_IVT_ID', 'Background', 1)
如果想设定物种信息,可用下面的方式设置,物种信息要用Taxonomy ID:
client.service.setCurrentSpecies('9606')
更多可用的函数见下图:
