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")+"\t"+listi[1]+"\t"+str(int(listi[1])+999999)+"\t"+str(float(listi[3])/10)+"\n")
   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")+"\t"+listx[1]+"\t"+str(int(listx[1])+499999)+"\t"+listx[3]+"\n")
      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))+"\t"+"\t".join(listlinks[0:3]).replace("chr","hs")+"\n")
            filesave.write("segdups"+str(listrandom.index(links))+"\t"+"\t".join(listlinks[5:8]).replace("chr","hs")+"\n")
         filesave.close()
      else:
         for links in listrandom:
            listlinks = links.split()
            filesave.write("segdups"+str(listrandom.index(links))+"\t"+"\t".join(listlinks[0:3])+"\n")
            filesave.write("segdups"+str(listrandom.index(links))+"\t"+"\t".join(listlinks[5:8])+"\n")
         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