我們很多Markers通過本地blast比對到多個物種的基因組上,來確定marker在不同物種基因組上的位置,為了方便查看結(jié)果,因此我寫了一個腳本,只需要輸入Markers文件和Blast結(jié)果所在的目錄,就能生成一個CSV結(jié)果文件,行名是Marker名稱,列名是物種名稱
代碼如下:
import numpy as np
import pandas as pd
import os
import csv
import argparse
import sys
###構(gòu)建參數(shù)傳遞
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")
###根據(jù)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)
###根據(jù)各個物種blast的結(jié)果提取物種名
species_list=[]
file_list=os.listdir(blast_result)
for file in file_list:
species=file.split(".")[0]###提出物種名
species_list.append(species)
###構(gòu)建輸出矩陣
data=np.zeros((len(marker_list),len(species_list)+1))
DF=pd.DataFrame(data)
###將所有的值都設(shè)置為NA
for col in range(0,len(species_list)+1):
for row in range(0,len(marker_list)):
DF.iloc[row,col]="NA"
##重置數(shù)據(jù)框的索引
rownames= {}
for i in range(0,len(marker_list)):
marker=marker_list[i]
rownames[i]=marker
DF.rename(index=rownames,inplace=True)##索引
###重置數(shù)據(jù)框的列名
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結(jié)果,修改矩陣中對應(yīng)的值
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文件