718 lines
26 KiB
Python
718 lines
26 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
# -*- coding: UTF-8 -*-
|
|||
|
|
|
|||
|
|
"""
|
|||
|
|
Created on: 2023-02-10
|
|||
|
|
@author: cjs
|
|||
|
|
# 用途：处理放化疗信息，重新构建脚本处理流程
|
|||
|
|
# 版本：1.0.0
|
|||
|
|
# 表格注意事项，杂合子顺序无影响，即AG和GA是一样的
|
|||
|
|
# 表格表头内容和顺序不能变更
|
|||
|
|
# 每个RS信息的第一个信息必须是野生型信息，未检测出的位点都按野生型计算
|
|||
|
|
# 最后编辑日期: 2023-02-10
|
|||
|
|
# 1.0.1 2023-03-10 药物只有一种型(LX、DF)的检测时，只判断一种型
|
|||
|
|
# 1.0.1 2023-03-10 药物的一个rs有多种型的检测，每种型都判读
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
from cjs_test.cjs_logger import Logger
|
|||
|
|
from openpyxl import load_workbook
|
|||
|
|
from glob import glob
|
|||
|
|
import datetime
|
|||
|
|
import traceback
|
|||
|
|
import os
|
|||
|
|
import sys
|
|||
|
|
import subprocess
|
|||
|
|
|
|||
|
|
# 全局参数
|
|||
|
|
Exe_Bin = ''
|
|||
|
|
Exe_Path = ''
|
|||
|
|
GLog = None
|
|||
|
|
StartTime = None
|
|||
|
|
# 命令参数
|
|||
|
|
Output_dir = r''
|
|||
|
|
Chemo_Path = ''
|
|||
|
|
Normal = ''
|
|||
|
|
Tumor = ''
|
|||
|
|
N_Bam = ''
|
|||
|
|
Bam_Path = ''
|
|||
|
|
Human_ref = '/dataseq/jmdna/database/genome/hg19/hg19.fa'
|
|||
|
|
Project = ''
|
|||
|
|
Pro_xlsx = ''
|
|||
|
|
# 运行变量
|
|||
|
|
Drug_Dd = {} # 记录rs对药物的影响
|
|||
|
|
Drug_Ls = [] # 继续处理药物的顺序
|
|||
|
|
Rs_Dd = {} # 构建bed文件所需的信息
|
|||
|
|
Rs_Gene = {} # rs对应的基因名称
|
|||
|
|
Error_Ls = ['Error_Type\tError_info\n']
|
|||
|
|
# 固定参数
|
|||
|
|
MIN_DEPTH = 30 # 检出的最低深度，低于这个不算检出
|
|||
|
|
MIN_RATE = 0.1 # 杂合的最低频率，低于这个补算杂合突变
|
|||
|
|
|
|||
|
|
RS_BED = 'rs.hg19.bed'
|
|||
|
|
# 表头信息可以左右平移，但是内容和顺序不能变动
|
|||
|
|
HEAD_S = '基因名称' # 需要处理列的开始
|
|||
|
|
HEAD_E = '证据等级' # 需要处理列的结束
|
|||
|
|
SNP_TYPES = ['LX', 'DF'] # 只解析的类型
|
|||
|
|
|
|||
|
|
# 可以运行的项目字典信息
|
|||
|
|
PROS_DD = {'160gene': '160.化疗位点注释及指导说明2023.02.02',
|
|||
|
|
'624gene': '624.化疗位点注释及指导说明2023.03.15.xlsx',
|
|||
|
|
'650gene': '650.panel化疗位点注释及指导说明.庞杰.20230403.xlsx'}
|
|||
|
|
|
|||
|
|
|
|||
|
|
def Exit_Print(pline=''):
|
|||
|
|
"""显示错误的信息，退出脚本."""
|
|||
|
|
print("Usage:")
|
|||
|
|
eline = "%s " % Exe_Bin
|
|||
|
|
eline += '-o Output_dir -p Project '
|
|||
|
|
eline += '[--n normal --t tumor --r Human_ref --b bam]'
|
|||
|
|
print(eline)
|
|||
|
|
if len(pline) > 0:
|
|||
|
|
print(pline)
|
|||
|
|
|
|||
|
|
if GLog:
|
|||
|
|
GLog.info('exit')
|
|||
|
|
GLog.close()
|
|||
|
|
sys.exit(0)
|
|||
|
|
|
|||
|
|
|
|||
|
|
# 处理运行参数
|
|||
|
|
def Get_Opts():
|
|||
|
|
"""获取运行的环境变量."""
|
|||
|
|
global Output_dir
|
|||
|
|
global Chemo_Path
|
|||
|
|
global Normal
|
|||
|
|
global Tumor
|
|||
|
|
global N_Bam
|
|||
|
|
global Bam_Path
|
|||
|
|
global Human_ref
|
|||
|
|
global Project
|
|||
|
|
global Pro_xlsx
|
|||
|
|
|
|||
|
|
global Exe_Bin
|
|||
|
|
global Exe_Path
|
|||
|
|
global GLog
|
|||
|
|
global StartTime
|
|||
|
|
|
|||
|
|
argvs = sys.argv
|
|||
|
|
file_real = os.path.realpath(argvs[0])
|
|||
|
|
Exe_Path = os.path.dirname(file_real)
|
|||
|
|
Exe_Bin = os.path.basename(file_real)
|
|||
|
|
StartTime = datetime.datetime.now()
|
|||
|
|
ymd = StartTime.__format__('%Y%m%d_%H%M%S')
|
|||
|
|
|
|||
|
|
pline = ''
|
|||
|
|
opt_normal = ''
|
|||
|
|
opt_tumor = ''
|
|||
|
|
opt_ref = ''
|
|||
|
|
N_Bam = ''
|
|||
|
|
if len(argvs) < 5:
|
|||
|
|
pline = '参数列表数量不对, %s' % argvs
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
|
|||
|
|
else:
|
|||
|
|
# 参数解析
|
|||
|
|
argv_index = 1
|
|||
|
|
for argv in argvs[1:]:
|
|||
|
|
argv_index += 1
|
|||
|
|
argv = argv.upper()
|
|||
|
|
|
|||
|
|
# 必选参数获取
|
|||
|
|
if argv == '-O':
|
|||
|
|
Output_dir = os.path.realpath(argvs[argv_index])
|
|||
|
|
# 新建日志
|
|||
|
|
if not os.path.exists(Output_dir):
|
|||
|
|
pline = '项目路径不存在: %s' % Output_dir
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
Chemo_Path = os.path.join(Output_dir, 'chemo')
|
|||
|
|
if not os.path.exists(Chemo_Path):
|
|||
|
|
os.makedirs(Chemo_Path)
|
|||
|
|
elif argv == '-P':
|
|||
|
|
Project = argvs[argv_index]
|
|||
|
|
if Project in PROS_DD:
|
|||
|
|
Pro_xlsx = os.path.join(Exe_Path, PROS_DD[Project])
|
|||
|
|
else:
|
|||
|
|
pline = '无法分析的项目类型:%s' % Project
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
|
|||
|
|
# 可选参数获得
|
|||
|
|
elif argv == '--N':
|
|||
|
|
opt_normal = argvs[argv_index]
|
|||
|
|
elif argv == '--T':
|
|||
|
|
opt_tumor = argvs[argv_index]
|
|||
|
|
elif argv == '--R':
|
|||
|
|
opt_ref = argvs[argv_index]
|
|||
|
|
elif argv == '--B':
|
|||
|
|
N_Bam = argvs[argv_index]
|
|||
|
|
|
|||
|
|
# 核对参数解析结果
|
|||
|
|
if Output_dir == '':
|
|||
|
|
pline = '未指定项目路径'
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
|
|||
|
|
# 解析noraml
|
|||
|
|
Bam_Path = os.path.join(Output_dir, 'alignment')
|
|||
|
|
bams = []
|
|||
|
|
if N_Bam == '':
|
|||
|
|
if opt_normal == '':
|
|||
|
|
bams = glob('%s/*rmdup.bam' % (Bam_Path))
|
|||
|
|
# 针对umi的项目，再寻找一次
|
|||
|
|
if len(bams) == 0:
|
|||
|
|
bams = glob('%s/*consensus.mapped.bam' % (Bam_Path))
|
|||
|
|
else:
|
|||
|
|
bams = glob('%s/*%s*rmdup.bam' % (Bam_Path, opt_normal))
|
|||
|
|
if len(bams) == 0:
|
|||
|
|
bams = glob('%s/*%s*consensus.mapped.bam' %
|
|||
|
|
(Bam_Path, opt_normal))
|
|||
|
|
|
|||
|
|
if len(bams) == 0:
|
|||
|
|
if opt_normal == '':
|
|||
|
|
pline = '项目路径未发现rmdup.bam文件'
|
|||
|
|
else:
|
|||
|
|
pline = '项目路径未发现名称为:%s的相关rmdup.bam文件' % opt_normal
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
else:
|
|||
|
|
bams = [N_Bam]
|
|||
|
|
|
|||
|
|
print('项目中的bams:')
|
|||
|
|
bam_size = 0
|
|||
|
|
for bam in bams:
|
|||
|
|
temp_base = os.path.basename(bam)
|
|||
|
|
temp_size = os.path.getsize(bam) / 1024 / 1024 / 1024
|
|||
|
|
print('bam_base:%s, bam_size:%0.2fGB' % (temp_base, temp_size))
|
|||
|
|
if N_Bam == '':
|
|||
|
|
N_Bam = bam
|
|||
|
|
bam_size = temp_size
|
|||
|
|
else:
|
|||
|
|
if temp_size < bam_size:
|
|||
|
|
bam_size = temp_size
|
|||
|
|
N_Bam = bam
|
|||
|
|
|
|||
|
|
# 选取文件大小最小的文件作为正常对照
|
|||
|
|
bam_base = os.path.basename(N_Bam)
|
|||
|
|
print('使用最小的bam文件:%s' % bam_base)
|
|||
|
|
bam_name = bam_base.split('.')[0]
|
|||
|
|
Normal = bam_name
|
|||
|
|
|
|||
|
|
# 输出文件tumor解析
|
|||
|
|
if opt_tumor == '':
|
|||
|
|
Tumor = Normal
|
|||
|
|
else:
|
|||
|
|
# 如果有normal，变更为normal名称
|
|||
|
|
# Tumor = opt_tumor
|
|||
|
|
if Normal != '':
|
|||
|
|
Tumor = Normal
|
|||
|
|
else:
|
|||
|
|
Tumor = opt_tumor
|
|||
|
|
|
|||
|
|
if opt_ref == '':
|
|||
|
|
print('参数中未解析出--R的参考基因组，使用默认文件:%s' % Human_ref)
|
|||
|
|
else:
|
|||
|
|
Human_ref = opt_ref
|
|||
|
|
|
|||
|
|
# 日志
|
|||
|
|
log_base = '%s_%s.log' % (Normal, ymd)
|
|||
|
|
log_full = os.path.join(Chemo_Path, log_base)
|
|||
|
|
GLog = Logger(log_full, mode='w')
|
|||
|
|
GLog.info('argvs, %s' % ' '.join(argvs))
|
|||
|
|
GLog.info('Normal:%s, Tumor:%s' % (Normal, Tumor))
|
|||
|
|
GLog.info('start')
|
|||
|
|
|
|||
|
|
|
|||
|
|
def Check_Rs_Type(r_drug, r_rs, r_type, d_type):
|
|||
|
|
"""检查杂合子，让杂合子对调后都有."""
|
|||
|
|
global Drug_Dd
|
|||
|
|
d_value = Drug_Dd[r_drug][r_rs][r_type][d_type]
|
|||
|
|
n_type = ''
|
|||
|
|
infos = []
|
|||
|
|
if r_type.find(r'/') > -1:
|
|||
|
|
infos = r_type.split(r'/')
|
|||
|
|
infos.insert(1, r'/')
|
|||
|
|
else:
|
|||
|
|
infos = list(r_type)
|
|||
|
|
for chr_str in infos:
|
|||
|
|
n_type = chr_str + n_type
|
|||
|
|
if n_type != r_type:
|
|||
|
|
if n_type not in Drug_Dd[r_drug][r_rs]:
|
|||
|
|
Drug_Dd[r_drug][r_rs][n_type] = {}
|
|||
|
|
# Drug_Dd[snp_drug][snp_rs].append(row_lines)
|
|||
|
|
Drug_Dd[r_drug][r_rs][n_type][d_type] = d_value
|
|||
|
|
|
|||
|
|
|
|||
|
|
def Check_Rs():
|
|||
|
|
"""检查项目所需的rs信息."""
|
|||
|
|
global Rs_Dd
|
|||
|
|
global Drug_Ls
|
|||
|
|
global Drug_Dd
|
|||
|
|
global Rs_Gene
|
|||
|
|
|
|||
|
|
bed_txt = os.path.join(Exe_Path, RS_BED)
|
|||
|
|
with open(bed_txt, 'r', encoding='utf8') as ff:
|
|||
|
|
for line in ff:
|
|||
|
|
line = line.replace('\r', '')
|
|||
|
|
line = line.replace('\n', '')
|
|||
|
|
if len(line) > 0 and not line.startswith('#'):
|
|||
|
|
lns = line.split('\t')
|
|||
|
|
rs_name = lns[3]
|
|||
|
|
if rs_name not in Rs_Dd:
|
|||
|
|
Rs_Dd[rs_name] = lns
|
|||
|
|
else:
|
|||
|
|
pline = '%s, rs有重复信息' % rs_name
|
|||
|
|
pline += '%s %s' % (Rs_Dd[rs_name], lns)
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
|
|||
|
|
wb = load_workbook(Pro_xlsx, read_only=True, data_only=False)
|
|||
|
|
wb_sheets = wb.sheetnames
|
|||
|
|
sheet = wb_sheets[0]
|
|||
|
|
ws1 = wb[sheet]
|
|||
|
|
sheet_rows = [row for row in ws1.rows] # 获取所有行
|
|||
|
|
# 表头处理
|
|||
|
|
head_row = sheet_rows[0]
|
|||
|
|
head_lines = []
|
|||
|
|
for cell in head_row:
|
|||
|
|
cell_str = str(cell.value)
|
|||
|
|
head_lines.append(cell_str)
|
|||
|
|
# print(head_lines)
|
|||
|
|
|
|||
|
|
head_start = head_lines.index(HEAD_S)
|
|||
|
|
head_end = head_lines.index(HEAD_E)
|
|||
|
|
end_pos = 0
|
|||
|
|
# 生成bed文件所需要的信息
|
|||
|
|
bed_rss = []
|
|||
|
|
bed_lines = []
|
|||
|
|
error_rss = []
|
|||
|
|
for row in sheet_rows[1:]:
|
|||
|
|
if end_pos == 1:
|
|||
|
|
break
|
|||
|
|
row_lines = []
|
|||
|
|
for col_index in range(head_start, head_end + 1):
|
|||
|
|
col_str = str(row[col_index].value)
|
|||
|
|
if col_str == 'None':
|
|||
|
|
end_pos = 1
|
|||
|
|
break
|
|||
|
|
row_lines.append(col_str)
|
|||
|
|
# 用药信息的字典构建
|
|||
|
|
snp_rs = row_lines[1]
|
|||
|
|
snp_drug = row_lines[5]
|
|||
|
|
snp_gene = row_lines[0]
|
|||
|
|
if snp_rs in Rs_Dd:
|
|||
|
|
Rs_Gene[snp_rs] = snp_gene
|
|||
|
|
|
|||
|
|
snp_infos = Rs_Dd[snp_rs]
|
|||
|
|
snp_chr = snp_infos[0]
|
|||
|
|
snp_start = snp_infos[1]
|
|||
|
|
snp_end = snp_infos[2]
|
|||
|
|
drug_type = row_lines[6].upper()
|
|||
|
|
# snp_type = row_lines[2].upper() # rs分型要不要统一大写，根据实际情况而定
|
|||
|
|
snp_type = row_lines[2]
|
|||
|
|
snp_score = int(row_lines[3])
|
|||
|
|
snp_res = row_lines[7]
|
|||
|
|
snp_level = row_lines[8]
|
|||
|
|
if snp_drug not in Drug_Dd:
|
|||
|
|
Drug_Dd[snp_drug] = {}
|
|||
|
|
Drug_Ls.append(snp_drug)
|
|||
|
|
if snp_rs not in Drug_Dd[snp_drug]:
|
|||
|
|
Drug_Dd[snp_drug][snp_rs] = {}
|
|||
|
|
|
|||
|
|
if snp_type not in Drug_Dd[snp_drug][snp_rs]:
|
|||
|
|
Drug_Dd[snp_drug][snp_rs][snp_type] = {}
|
|||
|
|
|
|||
|
|
# Drug_Dd[snp_drug][snp_rs].append(row_lines)
|
|||
|
|
Drug_Dd[snp_drug][snp_rs][snp_type][drug_type] = [snp_score,
|
|||
|
|
snp_res,
|
|||
|
|
snp_level]
|
|||
|
|
Check_Rs_Type(snp_drug, snp_rs, snp_type, drug_type)
|
|||
|
|
|
|||
|
|
if snp_rs not in bed_rss:
|
|||
|
|
bed_rss.append(snp_rs)
|
|||
|
|
snp_line = '\t'.join([snp_chr, snp_start, snp_end, snp_rs])
|
|||
|
|
bed_lines.append(snp_line + '\n')
|
|||
|
|
else:
|
|||
|
|
if snp_rs not in error_rss:
|
|||
|
|
error_rss.append(snp_rs)
|
|||
|
|
wb.close()
|
|||
|
|
|
|||
|
|
if len(error_rss) > 0:
|
|||
|
|
pline = '%s, 在%s文件中没有信息' % ('\t'.join(error_rss), RS_BED)
|
|||
|
|
Exit_Print(pline=pline)
|
|||
|
|
|
|||
|
|
# 写项目的bed文件
|
|||
|
|
bed_txt = '%s.bed' % Normal
|
|||
|
|
bed_full = os.path.join(Chemo_Path, bed_txt)
|
|||
|
|
with open(bed_full, 'w', encoding='utf-8') as ff:
|
|||
|
|
ff.writelines(bed_lines)
|
|||
|
|
return bed_full
|
|||
|
|
|
|||
|
|
|
|||
|
|
def Samtools_mpileup(df_bed):
|
|||
|
|
"""根据bed文件，分析文件中指定位置的突变信息,生成mpileup."""
|
|||
|
|
bam_mpileup = '%s.chemo.mpileup' % Normal
|
|||
|
|
bam_mpileup_full = os.path.join(Chemo_Path, bam_mpileup)
|
|||
|
|
sam_cmd = 'samtools mpileup -Bq 20 -Q 20 -a -f %s -l %s %s -o %s' % (
|
|||
|
|
Human_ref, df_bed, N_Bam, bam_mpileup_full)
|
|||
|
|
GLog.info('sam_cmd:%s' % sam_cmd)
|
|||
|
|
|
|||
|
|
print('Samtools_mpileup:%s' % bam_mpileup_full)
|
|||
|
|
ps = subprocess.Popen(sam_cmd,
|
|||
|
|
stdout=subprocess.PIPE,
|
|||
|
|
stderr=subprocess.STDOUT,
|
|||
|
|
shell=True)
|
|||
|
|
read_bytes = b''
|
|||
|
|
out_strs = ''
|
|||
|
|
while ps.poll() is None:
|
|||
|
|
out_std_b = ps.stdout.read(1)
|
|||
|
|
if out_std_b == b'\r' or out_std_b == b'\n':
|
|||
|
|
out_strs = read_bytes.decode().strip()
|
|||
|
|
if len(out_strs) > 0:
|
|||
|
|
GLog.info(out_strs)
|
|||
|
|
return 0
|
|||
|
|
|
|||
|
|
|
|||
|
|
def VarScan_mpileup2cns():
|
|||
|
|
"""mpileup-->VCF."""
|
|||
|
|
env_dist = os.environ
|
|||
|
|
# varscan 优先使用环境变量的版本
|
|||
|
|
varscan_exe = env_dist.get('VARSCAN')
|
|||
|
|
if varscan_exe == '':
|
|||
|
|
varscan_exe = '/dataseq/jmdna/software/VarScan.v2.4.2.jar'
|
|||
|
|
|
|||
|
|
bam_mpileup = '%s.chemo.mpileup' % Normal
|
|||
|
|
bam_mpileup_full = os.path.join(Chemo_Path, bam_mpileup)
|
|||
|
|
out_base = '%s.chemo.vcf' % Normal
|
|||
|
|
varscan_out = os.path.join(Chemo_Path, out_base)
|
|||
|
|
|
|||
|
|
varscan_cmd = 'java -jar %s ' % varscan_exe
|
|||
|
|
varscan_cmd += 'mpileup2cns %s ' % bam_mpileup_full
|
|||
|
|
varscan_cmd += '--min-var-freq 0.1 --min-freq-for-hom 0.9 '
|
|||
|
|
varscan_cmd += '--strand-filter 1 '
|
|||
|
|
varscan_cmd += '--output-vcf 1 >%s' % varscan_out
|
|||
|
|
GLog.info('varscan_cmd:%s' % varscan_cmd)
|
|||
|
|
print('VarScan_mpileup2cns:%s' % varscan_out)
|
|||
|
|
|
|||
|
|
ps = subprocess.Popen(varscan_cmd,
|
|||
|
|
stdout=subprocess.PIPE,
|
|||
|
|
stderr=subprocess.STDOUT,
|
|||
|
|
shell=True)
|
|||
|
|
read_bytes = b''
|
|||
|
|
out_strs = ''
|
|||
|
|
while ps.poll() is None:
|
|||
|
|
out_std_b = ps.stdout.read(1)
|
|||
|
|
if out_std_b == b'\r' or out_std_b == b'\n':
|
|||
|
|
out_strs = read_bytes.decode().strip()
|
|||
|
|
if len(out_strs) > 0:
|
|||
|
|
GLog.info(out_strs)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def Get_Drugs():
|
|||
|
|
"""处理VCF文件，核对结果，解析药物结果."""
|
|||
|
|
# 生成vcf字典
|
|||
|
|
vcf_dd = {}
|
|||
|
|
vcf_base = '%s.chemo.vcf' % Normal
|
|||
|
|
vcf_full = os.path.join(Chemo_Path, vcf_base)
|
|||
|
|
with open(vcf_full, 'r', encoding='utf-8') as ff:
|
|||
|
|
for line in ff:
|
|||
|
|
line = line.strip()
|
|||
|
|
if len(line) > 0 and not line.startswith('#'):
|
|||
|
|
lns = line.split('\t')
|
|||
|
|
snp_chr = lns[0]
|
|||
|
|
snp_end = lns[1]
|
|||
|
|
if snp_chr not in vcf_dd:
|
|||
|
|
vcf_dd[snp_chr] = {}
|
|||
|
|
vcf_dd[snp_chr][snp_end] = lns
|
|||
|
|
|
|||
|
|
# 处理药物的信息
|
|||
|
|
txt_head = '药物\t检测基因\t检测位点\t基因型\t证据等级\t用药提示\n'
|
|||
|
|
txt_res = [txt_head]
|
|||
|
|
drug_rd = {}
|
|||
|
|
err_rss = [] # 记录无法解析的rs信息
|
|||
|
|
for drug in Drug_Ls:
|
|||
|
|
for snp_rs in Drug_Dd[drug]:
|
|||
|
|
|
|||
|
|
snp_infos = Rs_Dd[snp_rs]
|
|||
|
|
snp_gene = Rs_Gene[snp_rs]
|
|||
|
|
snp_chr = snp_infos[0]
|
|||
|
|
snp_end = snp_infos[2]
|
|||
|
|
snp_ref = snp_infos[4]
|
|||
|
|
snp_alt = snp_infos[5]
|
|||
|
|
|
|||
|
|
# 解析VCF结果
|
|||
|
|
vcf_infos = vcf_dd.get(snp_chr).get(snp_end)
|
|||
|
|
vcf_mut_infos = vcf_infos[-1].split(':')
|
|||
|
|
vcf_alt = vcf_infos[3]
|
|||
|
|
vcf_mut = vcf_infos[4]
|
|||
|
|
vcf_rate = vcf_mut_infos[0]
|
|||
|
|
vcf_dep = 0
|
|||
|
|
vcf_fre = 0
|
|||
|
|
vcf_type = r''
|
|||
|
|
if len(vcf_mut_infos) > 10:
|
|||
|
|
vcf_dep = int(vcf_mut_infos[1])
|
|||
|
|
vcf_fre = float(vcf_mut_infos[6].replace(r'%', '')) / 100
|
|||
|
|
if vcf_dep < MIN_DEPTH:
|
|||
|
|
vcf_type = r'/'
|
|||
|
|
|
|||
|
|
# 判断VCF的检测结果
|
|||
|
|
mut_len = len(vcf_mut)
|
|||
|
|
alt_len = len(vcf_alt)
|
|||
|
|
if vcf_type != r'/': # 位点正常检出
|
|||
|
|
vcf_del = ''
|
|||
|
|
vcf_ins = ''
|
|||
|
|
if mut_len == 1 and alt_len == 1: # 单碱基处理
|
|||
|
|
if vcf_mut == '.': # 野生型
|
|||
|
|
vcf_type = '%s%s' % (vcf_alt, vcf_alt)
|
|||
|
|
else:
|
|||
|
|
if vcf_rate == r'1/1': # 纯合突变
|
|||
|
|
vcf_type = '%s%s' % (vcf_mut, vcf_mut)
|
|||
|
|
elif vcf_rate == r'0/1': # 杂合突变
|
|||
|
|
if vcf_fre < MIN_RATE: # 低频的杂合不算杂合
|
|||
|
|
vcf_type = '%s%s' % (vcf_alt, vcf_alt)
|
|||
|
|
else:
|
|||
|
|
vcf_type = '%s%s' % (vcf_alt, vcf_mut)
|
|||
|
|
snp_alt = snp_alt.split(',')[0] # 去除单碱基多种突变类型当作ins的可能
|
|||
|
|
# 针对del型的snp，野生型单独判断
|
|||
|
|
if len(snp_ref) > len(snp_alt):
|
|||
|
|
snp_wide = snp_ref.replace(snp_alt, '', 1)
|
|||
|
|
if len(snp_wide) == 1:
|
|||
|
|
vcf_type = '%s%s' % (snp_wide, snp_wide)
|
|||
|
|
else:
|
|||
|
|
vcf_type = '%s/%s' % (snp_wide, snp_wide)
|
|||
|
|
# 针对ins型的snp，野生型单独判断
|
|||
|
|
elif len(snp_ref) < len(snp_alt):
|
|||
|
|
vcf_type = '%s%s' % (snp_ref, snp_ref)
|
|||
|
|
else:
|
|||
|
|
# del型TC->T=C/DEL
|
|||
|
|
if alt_len > mut_len:
|
|||
|
|
vcf_del = vcf_alt.replace(vcf_mut, '', 1)
|
|||
|
|
if vcf_rate == r'1/1': # 纯合突变
|
|||
|
|
vcf_type = "del/del"
|
|||
|
|
elif vcf_rate == r'0/1': # 杂合突变
|
|||
|
|
if vcf_fre < MIN_RATE: # 低频的杂合不算杂合
|
|||
|
|
vcf_type = "%s/%s" % (snp_ref, snp_ref)
|
|||
|
|
else:
|
|||
|
|
vcf_type = vcf_del + r'/del'
|
|||
|
|
# Ins型
|
|||
|
|
elif alt_len < mut_len:
|
|||
|
|
vcf_ins = vcf_mut.replace(vcf_alt, '', 1)
|
|||
|
|
if vcf_rate == r'1/1': # 纯合突变
|
|||
|
|
vcf_type = "ins%s/ins%s" % (vcf_ins, vcf_ins)
|
|||
|
|
elif vcf_rate == r'0/1': # 杂合突变
|
|||
|
|
if vcf_fre < MIN_RATE: # 低频的杂合不算杂合
|
|||
|
|
vcf_type = "%s/%s" % (snp_ref, snp_ref)
|
|||
|
|
else:
|
|||
|
|
vcf_type = "%s/ins%s" % (vcf_alt, vcf_ins)
|
|||
|
|
GLog.info('RS:%s, VCF:%s, VCF_Fre:%s' % (
|
|||
|
|
snp_rs, vcf_type, vcf_fre))
|
|||
|
|
|
|||
|
|
# 个性化定制方案，针对多碱基突变的情况，每个rs号单独写
|
|||
|
|
if snp_rs == 'rs3064744': # 只考虑3种型，其它默认野生型
|
|||
|
|
# vcf_type:"C/insAT",vcf_ins:AT
|
|||
|
|
if len(vcf_ins) == 2: # 只考虑增加一个TA的情况
|
|||
|
|
if vcf_rate == r'0/1':
|
|||
|
|
if vcf_fre < MIN_RATE: # 低频的杂合不算杂合
|
|||
|
|
vcf_type = "TA[6]/TA[6]"
|
|||
|
|
else:
|
|||
|
|
vcf_type = "TA[6]/TA[7]"
|
|||
|
|
elif vcf_rate == r'1/1':
|
|||
|
|
vcf_type = "TA[7]/TA[7]"
|
|||
|
|
else:
|
|||
|
|
vcf_type = "TA[6]/TA[6]"
|
|||
|
|
else: # 野生型
|
|||
|
|
vcf_type = "TA[6]/TA[6]"
|
|||
|
|
GLog.info('RS:%s, VCF:%s, VCF_Fre:%s' % (
|
|||
|
|
snp_rs, vcf_type, vcf_fre))
|
|||
|
|
# 650基因
|
|||
|
|
elif snp_rs == 'rs45445694':
|
|||
|
|
# (CCGCGCCACTTGGCCTGCCTCCGTCCCG)3/(CCGCGCCACTTGGCCTGCCTCCGTCCCG)2
|
|||
|
|
if vcf_alt == 'A' and vcf_mut == '.': # 未突变
|
|||
|
|
vcf_alt = '(CCGCGCCACTTGGCCTGCCTCCGTCCCG)2'
|
|||
|
|
vcf_mut = '(CCGCGCCACTTGGCCTGCCTCCGTCCCG)2'
|
|||
|
|
# 缺失的强制变更为野生型
|
|||
|
|
elif vcf_alt == 'ACCGCGCCACTTGGCCTGCCTCCGTCCCG':
|
|||
|
|
if vcf_mut == 'A':
|
|||
|
|
vcf_alt = '(CCGCGCCACTTGGCCTGCCTCCGTCCCG)2'
|
|||
|
|
vcf_mut = '(CCGCGCCACTTGGCCTGCCTCCGTCCCG)2'
|
|||
|
|
vcf_type = '%s/%s' % (vcf_alt, vcf_mut)
|
|||
|
|
|
|||
|
|
# 判断有无检出
|
|||
|
|
snp_res = ''
|
|||
|
|
drug_type_key = ''
|
|||
|
|
snp_score = ''
|
|||
|
|
if vcf_type != r'/':
|
|||
|
|
# 获取用药指导
|
|||
|
|
if vcf_type not in Drug_Dd[drug][snp_rs]:
|
|||
|
|
if snp_rs not in err_rss:
|
|||
|
|
err_rss.append(snp_rs)
|
|||
|
|
eline = 'VCF结果无法解析\t%s,RS:%s' % (vcf_type, snp_rs)
|
|||
|
|
Error_Ls.append(eline + '\n')
|
|||
|
|
else:
|
|||
|
|
drug_type_dd = Drug_Dd[drug][snp_rs][vcf_type]
|
|||
|
|
drug_type_keys = list(drug_type_dd.keys())
|
|||
|
|
for drug_type_key in drug_type_keys:
|
|||
|
|
drug_infos = drug_type_dd[drug_type_key]
|
|||
|
|
snp_score = drug_infos[0]
|
|||
|
|
snp_res = drug_infos[1]
|
|||
|
|
snp_level = drug_infos[-1]
|
|||
|
|
snp_line = '\t'.join(
|
|||
|
|
[drug, snp_gene, snp_rs,
|
|||
|
|
vcf_type, snp_level, snp_res])
|
|||
|
|
txt_res.append(snp_line + '\n')
|
|||
|
|
|
|||
|
|
# 药物的LX和DF结果汇总
|
|||
|
|
if drug not in drug_rd:
|
|||
|
|
drug_rd[drug] = []
|
|||
|
|
# 收集药物的结果信息
|
|||
|
|
# '环磷酰胺': ['LX:1:疗效较好', 'LX:0:疗效居中', 'LX:-1:疗效较差']
|
|||
|
|
drug_value = '%s:%s:%s' % (drug_type_key,
|
|||
|
|
snp_score,
|
|||
|
|
snp_res)
|
|||
|
|
drug_rd[drug].append(drug_value)
|
|||
|
|
|
|||
|
|
else: # 未检出的也保留信息到结果文件，不是错误文件
|
|||
|
|
fs_type = list(Drug_Dd[drug][snp_rs].keys())[0]
|
|||
|
|
Drug_Fsd = Drug_Dd[drug][snp_rs][fs_type]
|
|||
|
|
for fs_key in Drug_Fsd:
|
|||
|
|
drug_infos = Drug_Fsd[fs_key]
|
|||
|
|
snp_level = drug_infos[2]
|
|||
|
|
snp_res = r'%s，无法判读' % fs_key
|
|||
|
|
snp_line = '\t'.join(
|
|||
|
|
[drug, snp_gene, snp_rs, vcf_type, snp_level, snp_res])
|
|||
|
|
txt_res.append(snp_line + '\n')
|
|||
|
|
# 药物的LX和DF结果汇总
|
|||
|
|
if drug not in drug_rd:
|
|||
|
|
drug_rd[drug] = []
|
|||
|
|
|
|||
|
|
if len(Error_Ls) > 1:
|
|||
|
|
# print(Error_Ls)
|
|||
|
|
error_base = '%s.drug.erros.txt' % Tumor
|
|||
|
|
error_full = os.path.join(Chemo_Path, error_base)
|
|||
|
|
with open(error_full, 'w', encoding='utf8') as ff:
|
|||
|
|
ff.writelines(Error_Ls)
|
|||
|
|
|
|||
|
|
# 生成文件
|
|||
|
|
drug_base = '%s.drug.infos.txt' % Tumor
|
|||
|
|
drug_full = os.path.join(Chemo_Path, drug_base)
|
|||
|
|
with open(drug_full, 'w', encoding='utf8') as ff:
|
|||
|
|
# 把相同rs的结果汇总，例如一个rs即知道疗效又指导毒副，那就写一行
|
|||
|
|
temp_drug = ''
|
|||
|
|
temp_gene = ''
|
|||
|
|
temp_rs = ''
|
|||
|
|
temp_line = ''
|
|||
|
|
line = ''
|
|||
|
|
for line in txt_res:
|
|||
|
|
lns = line.split('\t')
|
|||
|
|
line_drug = lns[0]
|
|||
|
|
line_gene = lns[1]
|
|||
|
|
line_rs = lns[2]
|
|||
|
|
if line_drug == temp_drug and \
|
|||
|
|
line_gene == temp_gene and line_rs == temp_rs:
|
|||
|
|
temp_line = temp_line.strip() + '，' + lns[-1]
|
|||
|
|
else:
|
|||
|
|
temp_drug = line_drug
|
|||
|
|
temp_gene = line_gene
|
|||
|
|
temp_rs = line_rs
|
|||
|
|
ff.write(temp_line)
|
|||
|
|
temp_line = line
|
|||
|
|
ff.write(temp_line)
|
|||
|
|
|
|||
|
|
# 解析是否推荐
|
|||
|
|
head = '药物名称\t疗效\t毒副\t推荐程度\t疗效和毒副总结\n'
|
|||
|
|
drug_lines = [head]
|
|||
|
|
drug_base = '%s.drug.res.txt' % Tumor
|
|||
|
|
drug_full = os.path.join(Chemo_Path, drug_base)
|
|||
|
|
# print(drug_rd['米托坦'])
|
|||
|
|
for drug in Drug_Ls:
|
|||
|
|
drug_lx = 0
|
|||
|
|
lx_find = 0
|
|||
|
|
drug_dx = 0
|
|||
|
|
dx_find = 0
|
|||
|
|
drug_res = ''
|
|||
|
|
drug_infos = ''
|
|||
|
|
snp_res_ls = drug_rd[drug]
|
|||
|
|
|
|||
|
|
if len(snp_res_ls) == 0:
|
|||
|
|
drug_line = '%s\t/\t/\t/\t/\n' % (drug)
|
|||
|
|
drug_lines.append(drug_line)
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
for snp_res in snp_res_ls:
|
|||
|
|
snp_ress = snp_res.split(':')
|
|||
|
|
# 检出
|
|||
|
|
if len(snp_ress) > 2:
|
|||
|
|
drug_type = snp_ress[0]
|
|||
|
|
drug_score = snp_ress[1]
|
|||
|
|
if drug_type.upper() == 'LX':
|
|||
|
|
drug_lx += int(drug_score)
|
|||
|
|
lx_find = 1
|
|||
|
|
elif drug_type.upper() == 'DF':
|
|||
|
|
drug_dx += int(drug_score)
|
|||
|
|
dx_find = 1
|
|||
|
|
|
|||
|
|
# 疗效和毒副汇总评估
|
|||
|
|
if drug_lx > 0:
|
|||
|
|
drug_infos = '疗效较好'
|
|||
|
|
elif drug_lx == 0 and lx_find == 1:
|
|||
|
|
drug_infos = '疗效一般'
|
|||
|
|
elif drug_lx < 0:
|
|||
|
|
drug_infos = '疗效较差'
|
|||
|
|
if drug_dx > 0:
|
|||
|
|
drug_infos += '，毒副较低'
|
|||
|
|
elif drug_dx == 0 and dx_find == 1:
|
|||
|
|
drug_infos += '，毒副一般'
|
|||
|
|
elif drug_dx < 0:
|
|||
|
|
drug_infos += '，毒副较高'
|
|||
|
|
if drug_infos.startswith('，'):
|
|||
|
|
drug_infos = drug_infos.replace('，', '', 1)
|
|||
|
|
|
|||
|
|
# print(drug, drug_lx, lx_find, drug_infos)
|
|||
|
|
# 最终用药判断
|
|||
|
|
if drug_lx > 0:
|
|||
|
|
if drug_dx >= 0:
|
|||
|
|
drug_res = '推荐'
|
|||
|
|
else:
|
|||
|
|
drug_res = '常规'
|
|||
|
|
elif drug_lx == 0:
|
|||
|
|
if drug_dx > 0:
|
|||
|
|
drug_res = '推荐'
|
|||
|
|
elif drug_dx == 0:
|
|||
|
|
drug_res = '常规'
|
|||
|
|
else:
|
|||
|
|
drug_res = '谨慎'
|
|||
|
|
else:
|
|||
|
|
if drug_dx > 0:
|
|||
|
|
drug_res = '常规'
|
|||
|
|
else:
|
|||
|
|
drug_res = '谨慎'
|
|||
|
|
# 个性化判断的药物再单独处理一下
|
|||
|
|
# 氟尿嘧啶、卡培他滨2种药疗效较好但毒副较高的结论为谨慎使用
|
|||
|
|
if drug_dx < 0:
|
|||
|
|
if drug == "氟尿嘧啶" or drug == "卡培他滨":
|
|||
|
|
drug_res = '谨慎'
|
|||
|
|
# 只做单项药物，注释掉没有做的项目，不再默认0，用斜号标识
|
|||
|
|
# if lx_find == 0:
|
|||
|
|
# drug_lx = r'/'
|
|||
|
|
# if dx_find == 0:
|
|||
|
|
# drug_dx = r'/'
|
|||
|
|
|
|||
|
|
drug_line = '%s\t%s\t%s\t%s\t%s\n' % (
|
|||
|
|
drug, drug_lx, drug_dx, drug_res, drug_infos)
|
|||
|
|
drug_lines.append(drug_line)
|
|||
|
|
with open(drug_full, 'w', encoding='utf8') as ff:
|
|||
|
|
ff.writelines(drug_lines)
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == '__main__':
|
|||
|
|
Get_Opts()
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
pass
|
|||
|
|
bed = Check_Rs()
|
|||
|
|
Samtools_mpileup(bed)
|
|||
|
|
VarScan_mpileup2cns()
|
|||
|
|
Get_Drugs()
|
|||
|
|
except BaseException:
|
|||
|
|
GLog.error(traceback.format_exc())
|
|||
|
|
print(traceback.format_exc())
|
|||
|
|
endtime = datetime.datetime.now()
|
|||
|
|
GLog.info('end')
|
|||
|
|
GLog.info('run time:%s seconds' % ((endtime - StartTime).seconds))
|
|||
|
|
GLog.close()
|