circos pipeline

# /usr/bin/env python
# coding=utf-8
###################################
# Author : yunkeli
# Version : 1.0(2015/6/20)
# E-mail : 1316014512@qq.com
###################################
import os
import argparse
import re
import random
def vcf_SNPdensity(snpvcffile,pathway):
print "this step is vcf to SNPdensity "
cmdvcftorate = "/home/liyunke/vcftools_0.1.12b/bin/vcftools --vcf "+snpvcffile+" --out " + pathway+"/SNPdensity100K --SNPdensity 1000000"
result_analysis=os.popen(cmdvcftorate)
print result_analysis.read()
def density(SNPdensityfile,pathway):
print "##############################"
print "this step is vcf to densitysplit cat "
fileopen = open(SNPdensityfile).readlines()[1:]
savename = pathway+"/"+"SNPdensity50K.snpden.new.txt"
filesave = open(savename,"w+")
for i in fileopen:
listi = i.split()
filesave.write(listi[0].replace("chr","hs")+" "+listi[1]+" "+str(int(listi[1])+999999)+" "+str(float(listi[3])/10)+" ")
filesave.close()
def densitysplit(SNPdensityfile,pathway):
print "##############################"
print "this step is vcf to densitysplit "
fileopen = open(SNPdensityfile).readlines()[1:]
namelist = []
for i in fileopen:
if i.split()[0] not in namelist:
namelist.append(i.split()[0])
for j in namelist:
savename = pathway+"/"+j.replace("chr","hs")+".snp.txt"
filesave = open(savename,"w+")
for x in fileopen:
listx = x.split()
if listx[0] == j:
filesave.write(listx[0].replace("chr","hs")+" "+listx[1]+" "+str(int(listx[1])+499999)+" "+listx[3]+" ")
filesave.close()
print "densitysplit ok"
def sv_split(svdensityfile,pathway):
print "##############################"
print "this step is vcf to sv_file split "
fileopen = open(svdensityfile).readlines()[1:]
namelist = []
for i in fileopen:
if i.split()[0] not in namelist:
namelist.append(i.split()[0])
for j in namelist:
listrandom = []
savename = pathway+"/"+j.replace("chr","hs")+".sv.txt"
filesave = open(savename,"w+")
for x in fileopen:
listx = x.split()
if listx[0] == j:
if listx[0] != listx[5]:
listrandom.append(x)
if len(listrandom) > 10:
slicelist = random.sample(listrandom, 10)
for links in slicelist:
listlinks = links.split()
filesave.write("segdups"+str(listrandom.index(links))+" "+" ".join(listlinks[0:3]).replace("chr","hs")+" ")
filesave.write("segdups"+str(listrandom.index(links))+" "+" ".join(listlinks[5:8]).replace("chr","hs")+" ")
filesave.close()
else:
for links in listrandom:
listlinks = links.split()
filesave.write("segdups"+str(listrandom.index(links))+" "+" ".join(listlinks[0:3])+" ")
filesave.write("segdups"+str(listrandom.index(links))+" "+" ".join(listlinks[5:8])+" ")
filesave.close()
def circos_config(npath,prefix):
oldconfig = "/home/liyunke/circos/sof/circos-0.67-7/work/pipeline/etc/config"
configopen = open(oldconfig).read()
f1 = re.sub("pathway",npath,configopen)
newconfig = "/home/liyunke/circos/sof/circos-0.67-7/work/pipeline/etc/"+prefix+".conf"
newconfigsave = open(newconfig,"w+")
newconfigsave.write(f1)
newconfigsave.close()
def main():
p = argparse.ArgumentParser(usage='./circos.pipline.py [--vcf] [--sv] [--prefix] [--outdir] ', description='circos snp sv')
p.add_argument('-v','--vcf', type=str, help='vcf file')
p.add_argument('-s','--sv', type=str, help='sv file')
p.add_argument('-p','--prefix', default="circostest",help='prefix or usrname')
p.add_argument('-o','--outdir', default="./", help='document directory')
args = p.parse_args()
prefix = args.prefix
vcffile = args.vcf
outdir = args.outdir
vcf_SNPdensity(vcffile,outdir)
SNPdensityfile = outdir+"/SNPdensity100K.snpden"
density(SNPdensityfile,outdir)
densitysplit(SNPdensityfile,outdir)
svdensityfile = args.sv
sv_split(svdensityfile,outdir)
circos_config(outdir,prefix)
cmdstr = "/home/liyunke/circos/sof/circos-0.67-7/bin/circos -conf /home/liyunke/circos/sof/circos-0.67-7/work/pipeline/etc/"+prefix+".conf --outputdir "+ outdir+" -outputfile "+prefix
result_analysis_circos =os.popen(cmdstr)
print result_analysis_circos.read()
rmcmd = "rm "+ outdir +"/hs*"
result_analysis_rm =os.popen(rmcmd)
print result_analysis_rm.read()
if __name__ == '__main__':
main()

原文地址:https://www.cnblogs.com/wipy/p/4600553.html