602panel_test/datafile_germline_v0_20230221_finish.py

import pandas as pd
import os,os.path
import xlrd
import re
import argparse

########germline针对的为白细胞对照
#保留条件如：频率>=30%、Func.refGene中intronic和ncRNA_exonic去掉、无突变注释信息的去掉、人群频率>=1%的去掉；
# 胚系注释中终止密码止为X表示，请跟体系保持一致，用*号。
# 胚系表头请按附近中示例调整一下顺序，VAF改为百分比的数据形式。
#去掉HLA相关基因突变
#输出按照基因名称进行排序

#for the end change
def rep(a):
	# print(a.group(0))
	st = re.sub("x","*",a.group(0),flags=re.I)
	st = re.sub("#","",st,flags=re.I)
	return st

#frameshift replace
def rep1(a):
	st = re.sub("fs\*.*", "", a.group(0), flags=re.I)
	st1 = re.split("(\d+)", st)
	if "," in a.group(0):
		st2 = st1[0] + st1[1] + "fs,"
	else:
		st2 = st1[0] + st1[1] + "fs"
	return st2

###当由多个基因构成时进行拆分
###当Gene.refGene由多个基因注释构成，且这几个基因是在我们的genelist中，我们需要的是将基因拆分，然后对应的AAchange信息保留
def gene_split(svdata):
	df_name = svdata['Gene.refGene'].str.split(';', expand=True)
	# 三、把行转列成列
	df_name = df_name.stack()
	# 四、重置索引，并删除多于的索引
	df_name = df_name.reset_index(level=1, drop=True)
	# 五、与原始数据合并
	df_name.name = 'df_name1'
	df_new = svdata.drop(['Gene.refGene'], axis=1).join(df_name)
	df_new.rename(columns={'df_name1': 'Gene.refGene'}, inplace=True)
	df_new.reset_index(drop=True, inplace=True)
	return df_new

###1.对主转录本数据的筛选。
#对'AAChange.refGene'的数据提取
#当'AAChange_select'为空时，
#1）当AAChange.refGene只有一个转录本的突变信息时，我们将这个转录本信息直接赋给AAChange_select
#2）当AAChange.refGene有多个转录本的突变信息时，选择这些转录本中最长的那个

#输入的是基因时
def transid_select_gene(gene):
	refseq_transdir = '/cgdata/liuxiangqiong/work62pancancer/Client/v0/script/refdata/hg19_refGene_transcript.bed'
	refseq = pd.read_table(refseq_transdir, sep='\t', usecols=[0, 1, 2, 6, 7], names=['chr', 'start', 'end', 'transid', 'gene'])
	refseq['trans_length'] = refseq['end'] - refseq['start']
	genetand = refseq[refseq['gene'] == gene]
	maxlength = genetand['trans_length'].max()
	trans_unique_infor = genetand[genetand['trans_length'] == maxlength]
	trans_unique_infor.reset_index(drop=True, inplace=True)
	trans_unique_id = trans_unique_infor.loc[0, 'transid']
	return trans_unique_id

#输入的时AAchange信息时.有一个或者多个转录本注释

def transid_select_AAchange(AAchange,AAchange_gene):
	AAchangedata0=AAchange.split(',')
	AAchangedata0_1=pd.DataFrame()
	for i in range(len(AAchangedata0)):
		AAchangedata0_1.loc[i,'AAchange']=AAchangedata0[i]
	AAchangedata1 =AAchangedata0_1['AAchange'].str.split(':', expand=True)
	AAchangedata1['transid'] = AAchangedata1[1].str.split('.', expand=True)[0]
	AAchangedata1.rename(columns={0: 'gene'}, inplace=True)
	AAchangedata2=AAchangedata1[AAchangedata1['gene']==AAchange_gene]
	refseq_transdir = '/cgdata/liuxiangqiong/work62pancancer/Client/v0/script/refdata/hg19_refGene_transcript.bed'
	refseq = pd.read_table(refseq_transdir, sep='\t', usecols=[0, 1, 2, 6, 7], names=['chr', 'start', 'end', 'transid', 'gene'])
	refseq['trans_length'] = refseq['end'] - refseq['start']
	##获得具有转录本编号的信息
	AAgene = pd.merge(AAchangedata2, refseq, on=['gene', 'transid'], how='inner')
	print(AAgene)
	#如果结果获得不为空则执行
	if len(AAgene) !=0:
		maxlength = AAgene['trans_length'].max()
		trans_unique_infor = AAgene[AAgene['trans_length'] == maxlength]
		trans_unique_infor.reset_index(drop=True, inplace=True)
		trans_unique_id = trans_unique_infor.loc[0, 'transid']
		geneinfor = AAchange_gene + ':' + trans_unique_id
		AAchange_select0 = [s for s in AAchangedata0 if geneinfor in s]
		AAchange_select=AAchange_select0[0]
	elif (len(AAgene) ==0) & (len(AAchangedata2)==1):
		AAchange_select =AAchange[0]
	elif (len(AAgene) ==0) & (len(AAchangedata2)>1):
		AAchange_select = AAchangedata0[0]
	return AAchange_select



####1.对AAchange_select的补充
def AAchangeselect_plus(germdata):
	for j in range(len(germdata)):
		#print(j)
		gene = germdata.loc[j, 'Gene.refGene']
		if pd.isnull(germdata.loc[j, 'AAChange_select']):
			germdata.loc[j, 'AAChange_select']=''
		if ( ',' in germdata.loc[j, 'AAChange_select']) and (pd.notnull(germdata.loc[j, 'AAChange_select'])):
			AAChange_selectinfor = germdata.loc[j, 'AAChange_select'].split(',')
			AAChange_selectinfor1 = [s for s in AAChange_selectinfor if gene+':' in s]
			germdata.loc[j, 'AAChange_select']=AAChange_selectinfor1
		AAchangedata = germdata.loc[j, 'AAChange.refGene']
		#print(AAchangedata)
		if (AAchangedata.count(',') == 0) & (AAchangedata != '.') & (germdata.loc[j, 'AAChange_select']==''):
			germdata.loc[j, 'AAChange_select'] = germdata.loc[j, 'AAChange.refGene']
		elif (AAchangedata.count(',') != 0) & (AAchangedata != '.') & (germdata.loc[j, 'AAChange_select']==''):
			gene_transid_unque = transid_select_AAchange(AAchangedata, gene)
			AAchange_infor = AAchangedata.split(',')
			#geneinfor = gene + ':' + gene_transid_unque
			geneinfor = gene_transid_unque
			AAchange_select = [s for s in AAchange_infor if geneinfor in s]
			if len(AAchange_select) != 0:
				germdata.loc[j, 'AAChange_select'] = AAchange_select[0]
			else:
				germdata.loc[j, 'AAChange_select'] = ''
	return germdata





#germdata1=AAchangeselect_plus(germdata)
		# 2.将c.注释的格式例如c.A8077T改为c.8077A>T
#germdata2=cdsmut_trans(germdata1)
#####2.将c.注释的格式例如c.A8077T改为c.8077A>T
#3.终止密码止原来用X表示，转换后用*表示
###4.对移码突变的转换,如由p.P1353Qfs*89改成p.P1353fs
def cdsmut_trans(germdata):
	newchange = germdata['AAChange_select'].str.split(":", expand=True)
	newchange['c1'] = newchange[3].str.split('.', expand=True)[1]
	for i in range(len(newchange)):
		if pd.notnull(newchange.loc[i, 'c1']):
			if ('ins' not in newchange.loc[i, 'c1']) & ('del' not in newchange.loc[i, 'c1']):
				new = 'c.' + re.findall("\d+\.?\d*", newchange.loc[i, 'c1'])[0] + re.sub("[0-9]+", '>',
																						 newchange.loc[i, 'c1'])
				newchange.loc[i, 'new_AA'] = new
			else:
				new = 'c.' + newchange.loc[i, 'c1']
				newchange.loc[i, 'new_AA'] = new
		else:
			newchange.loc[i, 'new_AA'] = ''
	germdata['new_select'] = newchange[0] + ':' + newchange[1] + ':' + newchange[2] + ':' + newchange['new_AA'] + ':' + newchange[4]
	# 终止密码止原来用X表示，转换后用*表示
	germdata['AAChange.refGene'] = (germdata['AAChange.refGene'] + "#").str.replace('p\..*?[,#]', rep)
	germdata['AAChange.refGene'] = germdata['AAChange.refGene'].str.replace('#$', "")
	germdata['new_select'] = (germdata['new_select'] + "#").str.replace('p\..*?[,#]', rep)
	germdata['new_select'] = (germdata['new_select']).str.replace('#$', "")
	# print(germline_filter4['AAChange_select'])
	###对移码突变的转换,如由p.P1353Qfs*89改成p.P1353fs
	germdata['AAChange.refGene'] = (germdata['AAChange.refGene'] + "#").str.replace('p\..*fs\*.*?[,#]', rep1)
	germdata['AAChange.refGene'] = germdata['AAChange.refGene'].str.replace('#$', "")
	germdata['new_select'] = (germdata['new_select'] + "#").str.replace('p\..*fs\*.*?[,#]', rep1)
	germdata['new_select'] = (germdata['new_select']).str.replace('#$', "")
	return germdata


####对AAchange列进一步,提取经过vep和annovar注释后的Otherinfo11这列信息
#并将氨基酸改变缩写改为单字母的简写，终止突变Ter改为*
def AAchange_new_infor(svdata):
	data_INFO_1 = svdata.loc[:, 'Otherinfo11'].str.split("|", expand=True)
	data_INFO_1.fillna("", inplace=True)
	result_num = 0
	for i in [i1 for i1 in range(1, data_INFO_1.shape[1], 41)]:
		data_INFO_2 = pd.DataFrame()
		data_INFO_2["AAChange.refGene_1"] = data_INFO_1.iloc[:, i + 2] + ":" + \
											data_INFO_1.iloc[:, i + 5] + ":" + \
											"exon" + data_INFO_1.iloc[:, i + 7].str.split("/", expand=True)[0] + ":" + \
											"intron" + data_INFO_1.iloc[:, i + 8].str.split("/", expand=True)[0] + ":"
		data_INFO_3 = data_INFO_1.iloc[:, i + 9].str.split(":", expand=True)
		data_INFO_4 = data_INFO_1.iloc[:, i + 10].str.split(":", expand=True)
		if data_INFO_3.shape[1] == 1:
			data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"] + data_INFO_3[0] + ":"
		else:
			data_INFO_3.fillna("", inplace=True)
			data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"] + data_INFO_3[1] + ":"
		if data_INFO_4.shape[1] == 1:
			data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"] + data_INFO_4[0]
		else:
			data_INFO_4.fillna("", inplace=True)
			data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"] + data_INFO_4[1]
		data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"].str.replace("exon:", "")
		data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"].str.replace("intron:", "")
		data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"].str.replace("::", ":")
		data_INFO_2["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"].str.replace(":$", "")
		if result_num == 0:
			svdata["AAChange.refGene_1"] = data_INFO_2["AAChange.refGene_1"] + ";"
			result_num = 1
		else:
			svdata["AAChange.refGene_1"] = svdata["AAChange.refGene_1"] + data_INFO_2["AAChange.refGene_1"] + ";"
	svdata["AAChange.refGene_1"] = svdata["AAChange.refGene_1"].str.replace(":;", "")
	svdata["AAChange.refGene_1"] = svdata["AAChange.refGene_1"].str.replace(";$", "")
	# 需要将氨基酸改变缩写改为单字母的简写，终止突变Ter改为*
	AA_path = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/AA_dict.txt'
	AAdata = pd.read_table(AA_path, sep='\t')
	for j in range(len(AAdata)):
		AAdata_abbre = AAdata.loc[j, 'Abbre_name']
		AAdata_short = AAdata.loc[j, 'Short_name']
		svdata['AAChange.refGene_1'] = svdata['AAChange.refGene_1'].str.replace(AAdata_abbre, AAdata_short)
	return svdata



##将vep和annovar注释结果进行整合
def vepannovar_merge(AAchange_raw,AAchange_new,AAchange_gene):
	'''
	vep and vardict anno merge
	输入数据均为字符串
	:param AAchange_raw: the AAchange.refgene from annovar anno of the data,one str。
	:param AAchange_new: the AAchange.refgene_1 from the vep anno of the data,one str
	:return: the modified AAchange,str
	'''
	###1.当AAchange_new中没有c.和p.信息时候，需要从AAchange_-raw中提取出来补全
	###1.1对AAchange_raw按照符合分列。
	AAchangeraw1 = AAchange_raw.split(',')
	###1.2对AAchane_new按照符号分列
	AAchangenew0 = AAchange_new.split(';')
	##从AAchane_new选择Gene.refGene在内的
	AAchangenew1 = [s for s in AAchangenew0 if AAchange_gene in s]
	if len(AAchangenew1) != 0:
		# 1.2.1循环分列后的AAchangenew1
		# 首先从割裂后的结果中，找到与AAchange_gene一致的突变
		# 然后从raw中查找对应的c.和p.信息。
		# 如果c.信息都不能找到，那么这个转录本对应的突变信息就设置为空。
		AAchangenew2 = ''
		for j in range(len(AAchangenew1)):
			#print(j)
			AAchangenew1_0 = AAchangenew1[j]
			#print(AAchangenew1_0)
			#如果信息中跨外显子的，那只提取外显子的
			if ('exon' in AAchangenew1_0) & ('intron' in AAchangenew1_0):
				AAchangenew1_1=re.sub('intron\d.*?:', "",AAchangenew1_0)
			else:
				AAchangenew1_1 =AAchangenew1_0
			#从AAchangeraw中把AAchangenew基因转录本对应的信息抓取出来
			id0 = AAchangenew1_1.split(':')
			id1 = id0[0] + ':' + id0[1].split('.')[0]
			select0 = [s for s in AAchangeraw1 if id1 in s]
			#print(select0)
			#3.当提取的信息不为空时，进行格式转换。如果抓取信息为空时，保留原来注释信息
			#c.注释的格式的转换：如c.A8077T改为c.8077A>T。由于AAchangeraw是annovar注释的，而AAchangenew时VEP注释
			if len(select0) != 0:
				# 如果信息中跨外显子的，那只提取外显子的
				if ('exon' in select0[0]) & ('intron' in select0[0]):
					select = re.sub('intron\d.*?:', "", select0[0])
				else:
					select = select0[0]
				select1 = select.split(':')
				#从分割后的list中体取出c.对应的信息
				cdsinfor0=[s for s in select1 if 'c.' in s]
				#如果能够找到c.信息，那么进行替换。如果找不到c.信息。那么新赋值的注释为空。我们只提取包含c.的，对于注释为n.的舍弃
				if len(cdsinfor0)!=0:
					changeinfro_select1=cdsinfor0[0].split(".")[1]
					if ('ins' not in changeinfro_select1) & ('del' not in changeinfro_select1)& ('dup' not in changeinfro_select1) & ('inv' not in changeinfro_select1):
						cdsinfor1 = 'c.' + re.findall("\d+\.?\d*", changeinfro_select1)[0] + re.sub("[0-9]+", '>',changeinfro_select1)
					else:
						cdsinfor1='c.'+changeinfro_select1
					# 从分隔后的list中提取楚p.对应的信息
					ppinfor0 = [s for s in select1 if 'p.' in s]
					if len(ppinfor0) != 0:
						ppinfor1 = ppinfor0[0]
					else:
						ppinfor1 = ''
					#获得新的注释
					annonew=id0[0] + ':' + id0[1]+':'+select1[2]+':'+cdsinfor1+':'+ppinfor1
				else:
					annonew =''
				# 4.对AAchangenew进行信息的填补
				# 如果AAchangenew没有c.信息时，那么新的赋值就用AAchangeraw里面
				if 'c.' not in AAchangenew1_1:
					AAchangenew2 = AAchangenew2 + annonew + ';'
				elif ('c.' in AAchangenew1_1) & ('p.' not in AAchangenew1_1):
					AAchangenew2 = AAchangenew2 + annonew + ';'
				else:
					AAchangenew2 = AAchangenew2 + AAchangenew1_1 + ';'
			else:
				AAchangenew2=AAchangenew2+AAchangenew1_1+';'
		AAchangenew3 = re.sub(r"(;)\1+", '', AAchangenew2)
		if AAchangenew3[-1] == ';':
			AAchangenew3_0 = AAchangenew3[:-1]
			if AAchangenew3_0[-1] ==':':
				AAchangenew4 = AAchangenew3_0[:-1]
			else:
				AAchangenew4 = AAchangenew3_0
		else:
			AAchangenew4 = AAchangenew3
	else:
		AAchangenew4 = AAchange_raw
	return AAchangenew4


###改，先对annovar注释的AAchange列进行格式转换
####对annovar注释的信息进行格式修改
def AAchange_annovar_transdata(AAchange_annovar):
	'''
	:param AAchange_annovar: 输入为annovar注释的AAchange.refgene列的一行
	:return: 经过校正后的
	'''
	if (AAchange_annovar!='.') & (AAchange_annovar!='UNKNOWN'):
		annovarsplit = AAchange_annovar.split(',')
		annovar_new1 = ''
		for j in range(len(annovarsplit)):
			label_cpoint = annovarsplit[j].find('c.')
			label_ppoint = annovarsplit[j].find('p.')
			if (label_cpoint != -1) & (label_ppoint != -1):
				changeinfro_select1 = annovarsplit[j][(label_cpoint + 2):(label_ppoint - 1)]
				aainfor = annovarsplit[j][label_ppoint:]
				if ('ins' not in changeinfro_select1) & ('del' not in changeinfro_select1) & (
						'dup' not in changeinfro_select1) & ('inv' not in changeinfro_select1):
					cdsinfor1 = 'c.' + re.findall("\d+\.?\d*", changeinfro_select1)[0] + re.sub("[0-9]+", '>',
																								changeinfro_select1)
				else:
					cdsinfor1 = 'c.' + changeinfro_select1
				# annovar_new00 = annovarsplit[j][:label_cpoint] + cdsinfor1 + ':' + aainfor
				# 对终止密码子的替换
				# annovar_new00 = annovar_new00 + "@"
				# annovar_new0 = re.sub("p\..+?[,@]", rep2, annovar_new00)
				if 'X' in aainfor:
					aainfor1 = aainfor.replace('X', '*')
				else:
					aainfor1 = aainfor
				annovar_new0 = annovarsplit[j][:label_cpoint] + cdsinfor1 + ':' + aainfor1
			elif (label_cpoint != -1) & (label_ppoint == -1):
				changeinfro_select1 = annovarsplit[j][(label_cpoint + 2):(label_ppoint - 1)]
				if ('ins' not in changeinfro_select1) & ('del' not in changeinfro_select1):
					cdsinfor1 = 'c.' + re.findall("\d+\.?\d*", changeinfro_select1)[0] + re.sub("[0-9]+", '>',
																								changeinfro_select1)
				else:
					cdsinfor1 = 'c.' + changeinfro_select1
				annovar_new0 = annovarsplit[j][:label_cpoint] + cdsinfor1
			else:
				annovar_new0 = annovarsplit[j]
			annovar_new1 = annovar_new0 + ',' + annovar_new1
		annovar_new_last = annovar_new1[:-1]
	else:
		annovar_new_last=AAchange_annovar
	return  annovar_new_last



def AAchange_merge(svdata):
	svdata.reset_index(drop=True, inplace=True)
	for i in range(len(svdata)):
		#print(i)
		svdata.loc[i, 'AAChange.refGene'] = AAchange_annovar_transdata(svdata.loc[i, 'AAChange.refGene'])
		AAchange_raw = svdata.loc[i, 'AAChange.refGene']
		AAchange_new = svdata.loc[i, 'AAChange.refGene_1']
		AAchange_gene = svdata.loc[i, 'Gene.refGene']
		exonfunc = svdata.loc[i, 'ExonicFunc.refGene']
		svdata.loc[i, 'AAChange.refGene_2'] = vepannovar_merge(AAchange_raw, AAchange_new,AAchange_gene)
		# 对AAChange.refGene_2进一步修改
		new2 = svdata.loc[i, 'AAChange.refGene_2']
		new21 = new2.split(';')
		if exonfunc == 'stopgain':
			# for stopgain,extract the mut including *
			new2_extract = [s for s in new21 if '*' in s]
			newchange = ';'.join(new2_extract)
		elif exonfunc == 'startloss':
			# for the startloss,extract the mut including ?
			new2_extract = [s for s in new21 if '?' in s]
			newchange = ';'.join(new2_extract)
		else:
			newchange = new2
		# 将%3D替换为=
		newchange_re = newchange.replace('%3D', '=')
		svdata.loc[i, 'AAChange.refGene_2'] = newchange_re
	return svdata


#######对转录本和突变表示形式的修改
####修改：去除化疗相关位点
###读入vep注释信息，方便对同义突变进一步筛选
def vepanno(inputpath,normalid):
	germdir = os.path.join(inputpath, '4Germline_unpair')
	# 读入vep注释
	vepdir = os.path.join(germdir, normalid + '.hg19_multianno.txt')
	if os.path.exists(vepdir) and os.path.getsize(vepdir) != 0:
		vepdata = pd.read_table(vepdir, sep='\t', header=0, low_memory=False)
		vepdata1 = AAchange_new_infor(vepdata)
		vepdata2 = AAchange_merge(vepdata1)
		#提取信息
		vepdata2['muttype'] = vepdata2['Otherinfo11'].str.split('|', expand=True)[1]
		vepdata3 = vepdata2[['Chr', 'Start', 'End', 'Ref', 'Alt', 'muttype', 'AAChange.refGene_1', 'AAChange.refGene_2']]
		vepdata3['Chr'] = vepdata3['Chr'].astype('str')
		vepdata3['Start'] = vepdata3['Start'].astype('str')
		vepdata3['End'] = vepdata3['End'].astype('str')
	return vepdata3

####改对unique的结果，根据vep的结果进一步校正。
####2.对于ExonicFunc.refGene为startloss选择结尾带？的突变信息；当为stopgain,选择下划线后面的信息。



#对AAChange_select进一步核对
def AAchange_select_new_test(germdata1):
	for i in range(len(germdata1)):
		aaannovar = germdata1.loc[i, 'AAChange_select'].split(':')
		aavep = germdata1.loc[i, 'AAChange.refGene_2'].split(';')
		if germdata1.loc[i,'ExonicFunc.refGene']=='startloss':
			result_aa = [s for s in aavep if '?' in s]
			germdata1.loc[i, 'AAChange_select'] = result_aa[0]
		elif germdata1.loc[i,'ExonicFunc.refGene']=='stopgain':
			result_aa = [s for s in aavep if '*' in s]
			germdata1.loc[i, 'AAChange_select'] = result_aa[0]
		else:
			seachgene = aaannovar[0] + ':' + aaannovar[1]
			result_aa = [s for s in aavep if seachgene in s]
			if len(result_aa) != 0:
				germdata1.loc[i, 'AAChange_select'] = result_aa[0]
	return germdata1


###对stopgain进一步修改
def AAchange_select_new_raw(germdata1):
	for i in range(len(germdata1)):
		print(i)
		aaannovar = germdata1.loc[i, 'AAChange_select'].split(':')
		if len(germdata1.loc[i, 'AAChange.refGene_2'])!=0:
			germdata1.loc[i, 'AAChange.refGene_2'] = germdata1.loc[i, 'AAChange.refGene_2'].replace(',', ';')
			aavep = germdata1.loc[i, 'AAChange.refGene_2'].split(';')
			if germdata1.loc[i, 'ExonicFunc.refGene'] == 'startloss':
				result_aa = [s for s in aavep if '?' in s]
				germdata1.loc[i, 'AAChange_select'] = result_aa[0]
			elif germdata1.loc[i, 'ExonicFunc.refGene'] == 'stopgain':
				result_aa = [s for s in aavep if '*' in s]
				# 对于特殊的修改:p.K173_Y174ins* -> Y174ins*
				if 'ins*' in result_aa[0]:
					splitinfor = result_aa.split(':')
					newaa0 = splitinfor[4]
					if '_' in newaa0:
						newaa1 = newaa0.split('_')
						numbefore = int(re.findall("\d+", newaa1[0])[0])
						numafter = int(re.findall("\d+", newaa1[1])[0])
						if numafter - numbefore == 1:
							Pnew = "p." + newaa1[1]
							unique_new = splitinfor[0] + ':' + splitinfor[1] + ':' + splitinfor[2] + ':' + splitinfor[
								3] + ':' + Pnew
							germdata1.loc[i, 'AAChange_select'] = unique_new
						else:
							germdata1.loc[i, 'AAChange_select'] = result_aa[0]
					else:
						germdata1.loc[i, 'AAChange_select'] = result_aa[0]
				else:
					germdata1.loc[i, 'AAChange_select'] = result_aa[0]
			else:
				seachgene = aaannovar[0] + ':' + aaannovar[1]
				result_aa = [s for s in aavep if seachgene in s]
				if len(result_aa) != 0:
					germdata1.loc[i, 'AAChange_select'] = result_aa[0]
	return germdata1

####改，当有多个的时候用唯一转录本

###对stopgain进一步修改
def AAchange_select_new(germdata1):
	for i in range(len(germdata1)):
		print(i)
		aaannovar = germdata1.loc[i, 'AAChange_select'].split(':')
		if len(germdata1.loc[i, 'AAChange.refGene_2'])!=0:
			germdata1.loc[i, 'AAChange.refGene_2'] = germdata1.loc[i, 'AAChange.refGene_2'].replace(',', ';')
			aavep = germdata1.loc[i, 'AAChange.refGene_2'].split(';')
			if germdata1.loc[i, 'ExonicFunc.refGene'] == 'startloss':
				result_aa = [s for s in aavep if '?' in s]
				germdata1.loc[i, 'AAChange_select'] = result_aa[0]
			elif germdata1.loc[i, 'ExonicFunc.refGene'] == 'stopgain':
				result_aa = [s for s in aavep if '*' in s]
				# 对于特殊的修改:p.K173_Y174ins* -> Y174ins*
				if 'ins*' in result_aa[0]:
					splitinfor = result_aa.split(':')
					newaa0 = splitinfor[4]
					if '_' in newaa0:
						newaa1 = newaa0.split('_')
						numbefore = int(re.findall("\d+", newaa1[0])[0])
						numafter = int(re.findall("\d+", newaa1[1])[0])
						if numafter - numbefore == 1:
							Pnew = "p." + newaa1[1]
							unique_new = splitinfor[0] + ':' + splitinfor[1] + ':' + splitinfor[2] + ':' + splitinfor[
								3] + ':' + Pnew
							germdata1.loc[i, 'AAChange_select'] = unique_new
						else:
							germdata1.loc[i, 'AAChange_select'] =(Transid_unique_search(','.join(result_aa))).replace(';','')
					else:
						germdata1.loc[i, 'AAChange_select'] = (Transid_unique_search(','.join(result_aa))).replace(';','')
				else:
					# 当不是p.K173_Y174ins*，从中选择唯一转录本所在的
					unqiue_trans_select =  (Transid_unique_search(','.join(result_aa))).replace(';','')
					if len(unqiue_trans_select)!=0:
						germdata1.loc[i, 'AAChange_select'] =unqiue_trans_select
					else:
						germdata1.loc[i, 'AAChange_select'] =result_aa[0]
			else:
				#首先选择唯一转录本
				unqiue_trans_select = (Transid_unique_search(','.join(aavep))).replace(';','')
				#如果没有获得唯一转录本的信息，采用AAchange_select原始的
				if len(unqiue_trans_select)!=0:
					germdata1.loc[i, 'AAChange_select'] =unqiue_trans_select
				else:
					seachgene = aaannovar[0] + ':' + aaannovar[1]
					result_aa = [s for s in aavep if seachgene in s]
					if len(result_aa) != 0:
						germdata1.loc[i, 'AAChange_select'] = result_aa[0]
			print(germdata1.loc[i, 'AAChange_select'] )
		#判断唯一转录本的区域是否和Func.refGene一致，如果不一致按照唯一转录本进行修改
		functype=germdata1.loc[i, 'Func.refGene'] 
		if 'intron' in germdata1.loc[i, 'AAChange_select']:
			func_unique='intronic'
		elif 'exon' in germdata1.loc[i, 'AAChange_select']:
			func_unique='exonic'
		else:
			func_unique='other'
		if (func_unique not in functype) & (func_unique!='other'):
			germdata1.loc[i, 'Func.refGene'] =func_unique
	return germdata1







###对唯一转录本的再次选择。
#首选在参考表中的转录本
###对于热门基因的转录本，采用主转录本
def Transid_unique_search(AAchange):
	transid_path = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/Select_RefSeq_HGNC_MANE.txt'
	transiddata = pd.read_table(transid_path, sep='\t')
	unique_trans_last = ''
	if AAchange != ' ':
		AAchange1 = AAchange.split(',')
		for j in range(len(AAchange1)):
			if AAchange1[j].count(":") == 4:
				AAchange1_0 = AAchange1[j].split(':')[1]
				AAchange1_exon = AAchange1[j].split(':')[2]
				AAchange1_cc = AAchange1[j].split(':')[3]
				AAchange1_pp = AAchange1[j].split(':')[4]
				transinfor = transiddata[(transiddata['RefSeq'] == AAchange1_0) | (transiddata['HGNC'] == AAchange1_0) | (transiddata['MANE'] == AAchange1_0)].reset_index(drop=True)
				if not transinfor.empty:
					unique_trans = transinfor.loc[0, 'Symbol'] + ':' + AAchange1_0 + ':' + AAchange1_exon + ':' + AAchange1_cc + ':' + AAchange1_pp+';'
				else:
					unique_trans = ''
				unique_trans_last = unique_trans_last + unique_trans
			elif AAchange1[j].count(":") == 3:
				AAchange1_0 = AAchange1[j].split(':')[1]
				AAchange1_exon = AAchange1[j].split(':')[2]
				AAchange1_cc = AAchange1[j].split(':')[3]
				transinfor = transiddata[
					(transiddata['RefSeq'] == AAchange1_0) | (transiddata['HGNC'] == AAchange1_0) | (
							transiddata['MANE'] == AAchange1_0)].reset_index(drop=True)
				if not transinfor.empty:
					unique_trans = transinfor.loc[
									   0, 'Symbol'] + ':' + AAchange1_0 + ':' + AAchange1_exon + ':' + AAchange1_cc+';'
				else:
					unique_trans = ''
				unique_trans_last = unique_trans_last + unique_trans
			elif AAchange1[j].count(":") == 2:
				AAchange1_0 = AAchange1[j].split(':')[1]
				AAchange1_exon = AAchange1[j].split(':')[2]
				transinfor = transiddata[
					(transiddata['RefSeq'] == AAchange1_0) | (transiddata['HGNC'] == AAchange1_0) | (
							transiddata['MANE'] == AAchange1_0)].reset_index(drop=True)
				if not transinfor.empty:
					unique_trans = transinfor.loc[0, 'Symbol'] + ':' + AAchange1_0 + ':' + AAchange1_exon+';'
				else:
					unique_trans = ''
				unique_trans_last = unique_trans_last + unique_trans
			else:
				unique_trans = ''
				unique_trans_last = unique_trans_last + unique_trans
	else:
		unique_trans_last = ''
	return unique_trans_last





###进行gnomAD_exome_EAS和STR区域的筛选
def snvfilter(inputpath,germline_filter1,sampleid):
	germdir = os.path.join(inputpath, '4Germline_unpair')
	STRdir = '/cgdata/liuxiangqiong/work62pancancer/Client/v0/script/refdata/UCSC_STR_LCR_simpleRepeat_genomicSuperDups.bed'
	germline_filter1.rename(columns={'AAChange_select': 'raw_AAChange_select'}, inplace=True)
	germline_filter1.rename(columns={'new_select': 'AAChange_select'}, inplace=True)
	germline_filter1.sort_values(by='Gene.refGene', inplace=True)
	# for the germline summary data
	germdata = '/cgdata/liuxiangqiong/work62pancancer/Client/v0/script/refdata/germline_counts_stastic20220815.txt'
	germref = pd.read_table(germdata, sep='\t', header=0)
	germref1 = germref[['mut_new', 'sample_counts']].drop_duplicates(keep='first')
	germref1.reset_index(drop=True, inplace=True)
	titlelist = germline_filter1.columns
	if len(germline_filter1) != 0:
		# select Func.refGene !='intronic'和'ncRNA_exonic'
		germline_filter2 = germline_filter1[
			(germline_filter1['Func.refGene'] != 'intronic') & (germline_filter1['Func.refGene'] != 'ncRNA_exonic')]
		if len(germline_filter2) != 0:
			# selct AAChange.refGene!='.'
			germline_filter3 = germline_filter2[germline_filter2['AAChange.refGene'] != '.']
			if len(germline_filter3) != 0:
				# select gnomAD_exome_EAS<=1%
				print('step6---select the vaf of genomAD')
				filterdata4_0 = germline_filter3[(germline_filter3['gnomAD_exome_EAS'] != '.')]
				# 6.1提取频率小于1%
				filterdata4_1 = filterdata4_0[(filterdata4_0['gnomAD_exome_EAS'].astype("float") <= 0.01)]
				# 6.2提取没有突变频率的
				filterdata4_2 = germline_filter3[(germline_filter3['gnomAD_exome_EAS'] == '.')]
				# 6.3将小于等于1%和没有突变频率的合并
				filterdata4_3 = filterdata4_2.append(filterdata4_1)
				germline_filter4 = filterdata4_3.reset_index(drop=True)
				germline_filter4.reset_index(drop=True, inplace=True)
				if len(germline_filter4) != 0:
					##remove the mutation in STR region
					germbed = germline_filter4[['Chr', 'Start', 'End', 'Gene.refGene']]
					outputname1 = os.path.join(germdir, sampleid + '_germline.bed')
					germbed.to_csv(outputname1, index=False, header=None, encoding='gbk', sep='\t')
					outputname2 = os.path.join(germdir, sampleid + '_germline_unSTR.bed')
					STRcmd = 'bedtools subtract -a ' + outputname1 + ' -b ' + STRdir + '>' + outputname2
					os.system(STRcmd)
					germline_unSTR = pd.read_table(outputname2, sep='\t', names=['Chr', 'Start', 'End', 'Gene.refGene'])
					if len(germline_unSTR) != 0:
						germline_unSTR[['Chr', 'Start', 'End']] = germline_unSTR[['Chr', 'Start', 'End']].astype('str')
						germline_filter4[['Chr', 'Start', 'End']] = germline_filter4[['Chr', 'Start', 'End']].astype(
							'str')
						germline_filter5 = pd.merge(germline_filter4, germline_unSTR,
													on=['Chr', 'Start', 'End', 'Gene.refGene'], how='inner')
						germline_filter5_1 = germline_filter5[titlelist]
						germline_filter5_1['mut_new'] = germline_filter5_1['Chr'] + '_' + germline_filter5_1[
							'Start'] + '_' + germline_filter5_1['End'] + '_' + germline_filter5_1['Ref'] + '_' + \
														germline_filter5_1['Alt']
						germline_filter_result = pd.merge(germline_filter5_1, germref1, on=['mut_new'], how='left')
						del germline_filter_result['mut_new']
						for i in range(len(germline_filter_result)):
							Func_refGene = germline_filter_result.loc[i, 'Func.refGene']
							ExonicFunc_refGene = germline_filter_result.loc[i, 'ExonicFunc.refGene']
							if (Func_refGene == 'splicing') & (ExonicFunc_refGene == '.'):
								germline_filter_result.loc[i, 'ExonicFunc.refGene'] = 'splicing'
						germline_filter_result.sort_values(by='Gene.refGene', inplace=True)
					else:
						print(sampleid + ' has no result when filter the STR')
						germline_filter_result = pd.DataFrame(columns=titlelist)
						germline_filter_result.loc[0, 'sampleid'] = sampleid
						germline_filter_result.loc[0, 'label'] = 'No result after filtering the STR'
					os.remove(outputname1)
					os.remove(outputname2)
				# print(germline_filter_result['AAChange_select'])
				else:
					print(sampleid + ' has no result when filter the genomeAD_EAS')
					germline_filter_result = pd.DataFrame(columns=titlelist)
					germline_filter_result.loc[0, 'sampleid'] = sampleid
					germline_filter_result.loc[0, 'label'] = 'No result after filtered by genomeAD_EAS'
		else:
			print(sampleid + ' has no result when filtering the intronic and ncRNA_exonic')
			germline_filter_result = pd.DataFrame(columns=titlelist)
			germline_filter_result.loc[0, 'sampleid'] = sampleid
			germline_filter_result.loc[0, 'label'] = 'No result after filtering the intronic and ncRNA_exonic'
	else:
		print(sampleid + ' has no result aftere vaf filtering')
		germline_filter_result = pd.DataFrame(columns=titlelist)
		germline_filter_result.loc[0, 'sampleid'] = sampleid
		germline_filter_result.loc[0, 'label'] = 'No result after filtering VAF>=0.3'
	return germline_filter_result



def germline_summary_control(inputpath,laneid,normalid,sampleid):
	germtitle_dir = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/germline_title.txt'
	#panel genelist
	panel_genelistdir = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/602gene_select.txt'
	panel_genelist = pd.read_table(panel_genelistdir, sep='\t', header=0, names=['Gene.refGene'])
	#druglist
	druglistdir = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/drug_genelist_20220829.txt'
	drug_genelist = pd.read_table(druglistdir, sep='\t', header=None, names=['Gene.refGene'])
	drug_genelist['drug_lable'] = 'Drug_gene'
	germtile = pd.read_table(germtitle_dir, sep='\t', header=0)
	titlelist_final = list(germtile.columns)
	titlelist_final.append('sample_counts')
	#对输出总表加上AAchange_1和AAchange_2
	#titlelist_summary=list(titlelist_final)
	#titlelist_summary.append('AAChange.refGene_1')
	#titlelist_summary.append('AAChange.refGene_2')
	germdir = os.path.join(inputpath, '4Germline_unpair')
	germdatafile=germdir+'/'+normalid+'.germ.xls'
	#读入vep注释
	vepdata1 = vepanno(inputpath,normalid)
	if os.path.exists(germdatafile) and os.path.getsize(germdatafile) !=0:
		#sampleid = normalid[:-2]
		# make the result dir
		result_dir = os.path.join(inputpath, 'resultfile')
		if not os.path.exists(result_dir):
			os.mkdir(result_dir)
		sample_dir = os.path.join(result_dir, sampleid)
		if not os.path.exists(sample_dir):
			os.mkdir(sample_dir)
		tempfile=os.path.join(inputpath,'tempfile')
		if not os.path.exists(tempfile):
			os.mkdir(tempfile)
		germ_tem_dir=os.path.join(tempfile,'germline')
		if not os.path.exists(germ_tem_dir):
			os.mkdir(germ_tem_dir)
		germdata_raw = pd.read_table(germdatafile, sep='\t')
		#先按照基因名排序
		germdata_raw.sort_values(by='Gene.refGene', inplace=True)
		germdata_raw.reset_index(drop=True, inplace=True)
		#修改列名
		germdata_raw.rename(columns={'CLNACC':'CLNALLELEID','CLNDBN':'CLNDN','CLNDSDB':'CLNDISDB','CLINVARREVSTATS':'CLNREVSTAT'},inplace=True)
		germdata_raw['Chr'] = germdata_raw['Chr'].astype('str')
		germdata_raw['Start'] = germdata_raw['Start'].astype('str')
		germdata_raw['End'] = germdata_raw['End'].astype('str')
		##对基因进行分列
		germdata01 = gene_split(germdata_raw)
		# 提取在gene列表中的基因
		germdata02 = pd.merge(germdata01, panel_genelist, on=['Gene.refGene'], how='inner')
		germdata02.insert(0, 'sampleid', sampleid)
		#去除化疗基因
		germdata03 = pd.merge(germdata02, drug_genelist, on=['Gene.refGene'], how='left')
		germdata = germdata03[germdata03['drug_lable'] != 'Drug_gene']
		germdata.reset_index(drop=True, inplace=True)
		###对CLINSIG进行分列
		if len(germdata[germdata['CLINSIG'].str.contains('\\[')]) != 0:
			germdata05 = germdata['CLINSIG'].str.split("[", expand=True)
			germdata05[1].fillna('.', inplace=True)
			germdata['CLNSIG'] = germdata05[0]
			germdata['CLNSIGCONF'] = germdata05[1].str.replace(']', '')
		else:
			germdata['CLNSIG'] = germdata['CLINSIG']
			germdata['CLNSIGCONF'] = '.'
		del germdata['CLINSIG']
		del germdata['drug_lable']
		print('clinsig')
		germdata.reset_index(drop=True, inplace=True)
		#1.add the AAchangedata
		germdata1=AAchangeselect_plus(germdata)
		# 2.将c.注释的格式例如c.A8077T改为c.8077A>T
		germdata2=cdsmut_trans(germdata1)
		# remove the HLA related gene
		germdata_result1 = germdata2[~germdata2['Gene.refGene'].str.contains('HLA')]
		germdata_result1.reset_index(drop=True, inplace=True)
		# select VAF>=0.3
		germline_filter1 = germdata_result1[germdata_result1['VAF'] >= 0.3]
		# VAF转为%
		germline_filter1['VAF'] = germline_filter1['VAF'].apply(lambda x: format(x, '.2%'))
		#进行gnomAD_exome_EAS和STR区域的筛选
		germline_filter_result0=snvfilter(inputpath,germline_filter1,sampleid)
		#对AAchange_select进行annovar和vep注释的合并进一步修改
		if len(vepdata1)!=0:
			germdata_vepanno0 = pd.merge(germline_filter_result0, vepdata1, on=['Chr', 'Start', 'End', 'Ref', 'Alt'], how='left')
			germdata_vepanno1 = germdata_vepanno0[~germdata_vepanno0['muttype'].str.contains('synonymous')]
			del germdata_vepanno1['muttype']
			germdata_vepanno1.reset_index(drop=True, inplace=True)
			#进行unique的筛选
			germdata_result=AAchange_select_new(germdata_vepanno1)
			#对移码突变的再次确认修改
			germdata_result['AAChange_select'] = (germdata_result['AAChange_select'] + "#").str.replace(
				'p\..*fs\*.*?[,#]', rep1)
			germdata_result['AAChange_select'] = germdata_result['AAChange_select'].str.replace('#$', "")
		else:
			germdata_result=germline_filter_result0
		#输出两种注释结果
		outputfile0=os.path.join(germ_tem_dir,laneid + '-' + sampleid + '.germline_raw_infor.txt')
		germdata_result.to_csv(outputfile0,sep='\t',index=False,header=True)
		#输出规定的列
		germdata_result.rename(columns={'AAChange.refGene':'AAChange.refGene_annovar','AAChange.refGene_2':'AAChange.refGene'},inplace=True)
		germline_filter_result=germdata_result[titlelist_final]
		#输出结果
		outputfile1 = os.path.join(sample_dir, laneid + '-' + sampleid + '.germline.xlsx')
		writer = pd.ExcelWriter(outputfile1)
		germline_filter_result.to_excel(writer, sheet_name='germline', index=False)
		writer.save()
		writer.close()
	else:
		print(normalid+' has no data or the result is null')
		germline_filter_result=pd.DataFrame()
		germline_filter_result.loc[0, 'sampleid'] = sampleid
		germline_filter_result.loc[0, 'label'] = 'No file'
		# make the temp dir
		temp_dir = os.path.join(inputpath, 'tempfile')
		if not os.path.exists(temp_dir):
			os.mkdir(temp_dir)
		bugfile_dir = os.path.join(temp_dir, 'bugfile')
		if not os.path.exists(bugfile_dir):
			os.mkdir(bugfile_dir)
		outputfile2 = os.path.join(bugfile_dir, laneid + '-' + sampleid + '.germline.nofile.log.txt')
		germline_filter_result.to_csv(outputfile2, index=False, header=True, encoding='gbk', sep='\t')


def germline_runmain(inputpath,normalid):
	laneid = inputpath.split('/')[-1].split('-')[-1]
	datasummarydir = os.path.join(inputpath, 'datasummary')
	isExists = os.path.exists(datasummarydir)
	if not isExists:
		os.makedirs(datasummarydir)
	sampledir = os.path.join(inputpath, laneid + '_sample_infor_label.txt')
	samplelist = pd.read_table(sampledir, sep='\t', header=0)
	sampledata = samplelist[samplelist['normal'] == normalid]
	sampledata.reset_index(drop=True, inplace=True)
	sampleid = sampledata.loc[0, 'samplename']
	germdir = os.path.join(inputpath, '4Germline_unpair')
	germdatafile = germdir + '/' + normalid + '.germ.xls'
	if os.path.exists(germdatafile) and os.path.getsize(germdatafile) != 0:
		germline_summary_control(inputpath, laneid, normalid, sampleid)
	else:
		print(normalid + ' is null')


if __name__=='__main__':
	parser = argparse.ArgumentParser(description='filter the germline')
	parser.add_argument('-i', '--inputpath', type=str, help='the path of lane')
	parser.add_argument('-s', '--normalid', type=str, help='the normal name of sample')
	args = parser.parse_args()
	Inputpath = args.inputpath
	Normalid = args.normalid
	germline_runmain(Inputpath,Normalid)



###for all sample in lane
def germlinesummary(inputpath,laneid):
	datasummarydir = os.path.join(inputpath, 'datasummary')
	isExists = os.path.exists(datasummarydir)
	if not isExists:
		os.makedirs(datasummarydir)
	sampledir = os.path.join(inputpath, laneid + '_sample_infor_label.txt')
	samplelist = pd.read_table(sampledir, sep='\t', header=0)
	germdatasummary = pd.DataFrame()
	for i in range(len(samplelist)):
		sampleid = samplelist.loc[i, 'samplename']
		normalid = samplelist.loc[i, 'normal']
		print(normalid)
		germdir = os.path.join(inputpath, '4Germline_unpair')
		germdatafile = germdir + '/' + normalid + '.germ.xls'
		if os.path.exists(germdatafile) and os.path.getsize(germdatafile) != 0:
			germre=germline_summary_control(inputpath, laneid, normalid,sampleid)
			germdatasummary =germdatasummary.append(germre)
		else:
			print(normalid+' is null')
			continue
	outputname = datasummarydir + '/' + laneid + '_table6_germlinne_datasummary.txt'
	germdatasummary.to_csv(outputname, index=False, header=True, encoding='gbk', sep='\t')