#! /usr/bin/env python3 import argparse import logging import os import pandas as pd import pysam class ChemoRun: def __init__(self, database, probe, cancer, project, output_dir, name, vcf): self.database = database self.probe = probe self.output_dir = output_dir self.name = name self.vcf = vcf self.cancer = cancer self.project = project @staticmethod def read_info(project): info_path = os.path.join(os.environ.get('DATABASE'), 'info.csv') read_res = pd.read_csv(info_path) read_res = read_res.fillna('.') read_res = read_res[read_res['project'] == project] gene_list = list() if any(read_res['chemotherapy_drug'].values): genecontent = read_res['chemotherapy_drug'].values[0] if genecontent != '.': gene_list = genecontent.split('/') if not gene_list: raise UserWarning(f'化疗检测基因为空！{project} ,请查看参数project是否正确！') return gene_list @staticmethod def get_drug_plan(x, drugsum): tlist = x.split('+') tdeslist = list() for tdes in tlist: if tdes.strip() in drugsum: t1_des = drugsum[tdes.strip()] tdeslist.append(t1_des) if '慎用' in tdeslist or '谨慎' in tdeslist: return '慎用' elif '推荐' in tdeslist: return '推荐' elif '常规' in tdeslist: return '可选' else: return '可选' @staticmethod def build_record_dict(vcf_file): vcfdata = pysam.VariantFile(vcf_file) records = {} for record in vcfdata: chrom = record.chrom pos = record.pos records[(chrom, pos)] = record return records def parsedata(self): """ 处理chemo_drug_rsid_snppos.xlsx, 并生成对应bed文件 """ data = pd.read_excel(self.database, None) gene_list = self.read_info(self.project) drug_rsid = data['drug_rsid'] drug_combine = data['drug_combine'] drug_rsid.fillna('.', inplace=True) drug_combine.fillna('.', inplace=True) drug_rsid_data = drug_rsid[ (drug_rsid['source'].str.contains(self.probe) & (drug_rsid['genename'].apply(lambda x: x in gene_list)))] drug_combine_data = drug_combine[drug_combine['treatResult'].apply(lambda x: x in self.cancer)] if drug_rsid_data.empty: logging.error(f'无此项目！{self.probe} ,请查看参数probe是否正确！') raise UserWarning(f'无此项目！{self.probe} ,请查看参数probe是否正确！') # 生成bed文件 beddata = drug_rsid_data[['chr', 'start', 'end', 'rsid']].drop_duplicates() beddata.rename(columns={'chr': '#chr'}, inplace=True) bedpath = os.path.join(self.output_dir, f'{self.probe}_chemo.bed') beddata.to_csv(bedpath, sep='\t', index=False) return drug_rsid_data, drug_combine_data, bedpath def parse_vcf(self, vcfpath, drug_rsid_data, drug_combine_data): records = self.build_record_dict(vcfpath) beddata = drug_rsid_data[['chr', 'start', 'end', 'rsid']].drop_duplicates() fliterdata = pd.DataFrame() for _, row in beddata.iterrows(): chrom = row['chr'] end = row['end'] fliter = pd.DataFrame() if (chrom, end) not in records: fliter = pd.concat([fliter, drug_rsid_data[ (drug_rsid_data['chr'] == chrom) & (drug_rsid_data['end'] == end) & (drug_rsid_data['genotype'] == '0/0') ]]) else: record = records[(chrom, end)] ref = record.ref alt = record.alts[0] gt = '/'.join(list(map(str, sorted(record.samples.get(record.samples.keys()[0]).get('GT'))))) fliter = pd.concat([fliter, drug_rsid_data[ (drug_rsid_data['chr'] == chrom) & (drug_rsid_data['end'] == end) & (drug_rsid_data['ref'] == ref) & (drug_rsid_data['alt'] == alt) & (drug_rsid_data['genotype'] == gt) ]]) if fliter.empty: raise UserWarning( 'chr: %s , end: %s 数据库未能匹配, 野生型0/0也未能匹配' % (chrom, end)) fliterdata = pd.concat([fliterdata, fliter]) # 生成过滤之后文件 respath = os.path.join(self.output_dir, f'{self.name}.chemo.res.csv') if not fliterdata.empty: fliterdata['probe'] = self.probe fliterdata.to_csv(respath, sep='\t', index=False) # 分类汇总 同位点，药物合并 drug.infos.txt drugrsid = fliterdata[['drugname', 'genename', 'rsid', 'result', 'level', 'tips', 'drugsort']] drugrsid = drugrsid.drop_duplicates() resdrugrsid = drugrsid.groupby(['drugname', 'genename', 'rsid', 'result', 'level', 'drugsort'])['tips'].agg( ','.join).reset_index() resdrugrsid.rename(columns= {'drugname': '药物', 'genename': '检测基因', 'rsid': '检测位点', 'result': '基因型', 'level': '证据等级', 'tips': '用药提示'}, inplace=True) resdrugrsid = resdrugrsid.sort_values(by=['drugsort', '药物', '检测基因']) resdrugrsid.to_csv(os.path.join(self.output_dir, f'{self.name}.drug.infos.txt'), index=False, sep='\t') # 药物 药物疗效 推荐程度合并 drug.res.txt drugtypesum = fliterdata[['drugname', 'drugtype', 'rsid', 'weights']] drugtypesum = drugtypesum.drop_duplicates() drugtyperes = list() drugsum = dict() for drug, drugdata in drugtypesum.groupby('drugname'): tipsnum = drugdata.groupby(['drugtype']).agg({'weights': 'sum'}).to_dict('index') sumlist = list() if 'LX' in tipsnum: LX = tipsnum['LX']['weights'] if LX > 0: lxdes = '疗效较好' lxnum = 1 elif LX == 0: lxdes = '疗效一般' lxnum = 0 else: lxdes = '疗效较差' lxnum = -1 sumlist.append(lxdes) else: LX = 0 lxnum = 0 if 'DF' in tipsnum: DF = tipsnum['DF']['weights'] if DF > 0: dfdes = '毒副较低' dfnum = 1 elif DF == 0: dfdes = '毒副一般' dfnum = 0 else: dfdes = '毒副较高' dfnum = -1 sumlist.append(dfdes) else: DF = 0 dfnum = 0 # 评价方式 疗效 1 0 -1, 毒副 1 0 -1 ，可形成9宫格 sumnum = lxnum + dfnum if sumnum > 0: sumdes = '推荐' elif sumnum == 0: sumdes = '常规' else: sumdes = '谨慎' # 特别药物处理 if (drug == "氟尿嘧啶" or drug == "卡培他滨") and DF < 0: sumdes = '谨慎' drugtyperes.append(dict( 药物名称=drug, 疗效=LX, 毒副=DF, 推荐程度=sumdes, 疗效和毒副总结=','.join(sumlist) )) drugsum[drug] = sumdes # 报告中展示药物有顺序 drugsort = fliterdata[['drugname', 'drugsort']].drop_duplicates() drugsort_dict = drugsort.set_index('drugname')['drugsort'].to_dict() drugtyperes_sort = sorted(drugtyperes, key=lambda x: ( drugsort_dict[x['药物名称']] if x['药物名称'] in drugsort_dict else 100, x['药物名称'])) # 生成数据 pd.DataFrame(drugtyperes_sort).to_csv(os.path.join(self.output_dir, f'{self.name}.drug.res.txt'), index=False, sep='\t') # 联合用药 if not drug_combine_data.empty: drug_combine_data['临床提示'] = drug_combine_data['用药方案'].apply(self.get_drug_plan, args=(drugsum,)) drug_combine_data.to_csv(os.path.join(self.output_dir, f'{self.name}.chemo.comb.txt'), index=False, sep='\t') def run(self): try: drug_rsid_data, drug_combine_data, bedpath = self.parsedata() self.parse_vcf(self.vcf, drug_rsid_data, drug_combine_data) except UserWarning as e: raise UserWarning(e) if __name__ == "__main__": parser = argparse.ArgumentParser(description="Chemotherapy Process Script") parser.add_argument('-d', '--database', help="Path to chemo_drug's database", required=True) parser.add_argument('-probe', '--probe', help="Probe name", required=True) parser.add_argument('-n', '--name', help="Name for sample", required=True) parser.add_argument('-v', '--vcf', help="germline vcf", required=True) parser.add_argument('-c', '--cancer', help="cancer", required=True) parser.add_argument('-p', '--project', help="Project", required=True) parser.add_argument('-o', '--output_dir', help="Output directory, default ./", default='') args = parser.parse_args() chemorun = ChemoRun(args.database, args.probe, args.cancer, args.project, args.output_dir, args.name, args.vcf) chemorun.run()