將Marker與多個物種blast結(jié)果輸出位CSV文件的python腳本

我們很多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文件
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容