snp关联分析及提取相关基因

aat.py

#coding:utf-8

import sys

#处理permutaion的阈值 返回SNPlist
def promusnp(mperm,p=0.001):
    f=open(mperm,"r")
    res=[]
    total=[]
    line=f.readline()
    for line in f.readlines():
        cons=line.split()
        if float(cons[2])<p:
            res.append(cons[1])
    f.close()
    return res

#将SNP的信息写入文件
def assoinfo(assoc,snps,write):
    f=open(assoc,"r")
    wrlins=[]
    line=f.readline()
    for line in f.readlines():
        cons=line.split()
        if cons[1] in snps:
            wrlins.append("\t".join([cons[1],cons[2],cons[8],cons[9],cons[11],cons[12]])+"\n")
    f.close()
    f=open(write,"w")
    f.writelines(wrlins)
    f.close()

#从SNP得基因ID
def getids(genefile,snps):
    f=open(genefile,"r")
    res=[]
    line=f.readline()
    for line in f.readlines():
        cons=line.split('\t')
        if cons[1] in snps:
            if len(cons)<5:continue
            for atom in cons[4:]:
                res.append(atom)
    f.close()
    return set(res)

#从基因ID得名称
def getgenes(info,ids,write):
    f=open(info,"r")
    res=[]
    for line in f.readlines():
        line=line.replace("\r",'')
        cons=line.split('\t')
        if cons[0] in ids:
            res.append(cons[1])
    f.close()
    f=open(write,"w")
    f.writelines(res)
    f.close()
    
snps=promusnp("chr6-1soga.assoc.mperm")
assoinfo("chr6-1soga.assoc",snps,"snpinfos")
geneids=getids("chr6_SNP2Gene.txt",snps)
getgenes("human_gene_info.txt",geneids,"genelist")