我们很多Markers通过本地blast比对到多个物种的基因组上,来确定marker在不同物种基因组上的位置,为了方便查看结果,因此我写了一个脚本,只需要输入Markers文件和Blast结果所在的目录,就能生成一个CSV结果文件,行名是Marker名称,列名是物种名称
代码如下:
import numpy as np
import pandas as pd
import os
import csv
import argparse
import sys
###构建参数传递
pars=argparse.ArgumentParser()
pars.add_argument("-marker_file",required=True,help="The Markerfile which formart is Fasta")
pars.add_argument("-blast_result",required=True,help="The Folders of blast results for each species")
args=pars.parse_args()
if args.marker_file:
marker_index=sys.argv.index("-marker_file")
marker_file=sys.argv[marker_index+1]
else:
raise Exception("Please input the markers file")
if args.blast_result:
result_index=sys.argv.index("-blast_result")
blast_result=sys.argv[result_index+1]
else:
raise Exception("Please input The Folders of blast results for each species")
###根据Makers文件提取Marker name并将Marker序列储存
marker_list=[]
seq_list=[]
with open(marker_file) as f:
for line in f:
if line.startswith(">"):
marker=line.replace(">","").split()[0]
marker_list.append(marker)
else:
seq=line.replace("\n","").strip()###删除空白字符
seq_list.append(seq)
###根据各个物种blast的结果提取物种名
species_list=[]
file_list=os.listdir(blast_result)
for file in file_list:
species=file.split(".")[0]###提出物种名
species_list.append(species)
###构建输出矩阵
data=np.zeros((len(marker_list),len(species_list)+1))
DF=pd.DataFrame(data)
###将所有的值都设置为NA
for col in range(0,len(species_list)+1):
for row in range(0,len(marker_list)):
DF.iloc[row,col]="NA"
##重置数据框的索引
rownames= {}
for i in range(0,len(marker_list)):
marker=marker_list[i]
rownames[i]=marker
DF.rename(index=rownames,inplace=True)##索引
###重置数据框的列名
columns={}
columns[0]='Sequences'
for i in range(1,len(species_list)+1):
species=species_list[i-1]
columns[i]=species
DF.rename(columns=columns,inplace=True)##列名
###将marker的序列输入
DF.loc[:,"Sequences"]=seq_list
###读取每个物种的blast结果,修改矩阵中对应的值
for file in file_list:
species = file.split(".") [0] ###提出物种名
with open(r"/projects01/DS20082800001/05.markerPosition/uniq_result/%s" % file) as f:
for line in f:
marker=line.split("\t")[0]
chr=line.split("\t")[1]
pos1=int(line.split("\t")[8])
pos2=int(line.split("\t")[9])
if(pos1<pos2):
start=str(pos1)
end=str(pos2)
else:
start=str(pos2)
end=str(pos1)
value=chr+":"+start+"-"+end
DF.loc[marker,species]=value
pd.DataFrame.to_csv(DF,"result.csv")###将矩阵输出为csv文件