602panel_test/datafile_QC_v0_20230406_finish.py

import sys,collections,math,os,os.path,re
from collections import Counter
import pandas as pd
from pandas.core.frame import DataFrame
import argparse
import xlrd
import subprocess
from time import time
from tqdm import tqdm
import threading
from multiprocessing.pool import ThreadPool #线程
from openpyxl import *
from openpyxl.styles import PatternFill, Border, Side, Alignment, Protection, Font






####excel颜色标签
# file:文件路径或workbook对象，第一次运行函数后会返回一个workbook对象，后续可直接调用这个返回的对象进行其他操作
# column_1:要操作的excel列,A-Z，比如第1列在excel中为A，第5列是E
# type_1:对比类型，1是>,2是<,3是==
# conditions_1:筛选条件,如type_1=1,conditions_1=10，则对大于等于10的值进行操作
# color_1:设置符合条件的值得颜色，可填red，blue，green，black四种
# reverse=False:设置不符合条件的值得颜色，默认是黑色，可填red，blue，green，black四种
# head=True:表中是否包含列名，默认为True
# percentage=False:是否是百分比数值，默认为False，当数值是百分比时设置为True
def excel_font_color(file, column_1, type_1, conditions_1,color_1, reverse=False, head=True, percentage=False):
	if type(file) == workbook.workbook.Workbook:
		wb2 = file
		ws2 = wb2.active
	else:
		wb2 = load_workbook(file)
		ws2 = wb2.active
	color = {"red": "ff000c", "green": "04ff00", "blue": "0006ff", "black":"000000"}
	bold_itatic_24_font_2 = Font(color=color[color_1])
	bold_itatic_24_font_3 = False
	if reverse:
		bold_itatic_24_font_3 = Font(color=color[reverse])
	num = 0
	if head:
		num = 1
	for i in ws2[column_1]:
		if num == 1:
			num = 0
			continue
		if num == 0:
			if "%" in str(i.value):
				print("{}列含有%,注意设置percentage=True,如不需设置则忽略此警告".format(column_1))
				num = 2
		# 带有百分比的值读取时值类型可能存在不同，需要分开设置
		if percentage:
			if i.number_format == "General":
				vAlue = float(i.value[:-1])/100
			else:
				vAlue = i.value
		else:
			vAlue = i.value
		# 数据加颜色
		if type_1 == 1:
			if vAlue > conditions_1:
				i.font = bold_itatic_24_font_2
			else:
				if bold_itatic_24_font_3:
					i.font = bold_itatic_24_font_3
		if type_1 == 2:
			if vAlue < conditions_1:
				i.font = bold_itatic_24_font_2
			else:
				if bold_itatic_24_font_3:
					i.font = bold_itatic_24_font_3
		if type_1 == 3:
			if vAlue == conditions_1:
				i.font = bold_itatic_24_font_2
			else:
				if bold_itatic_24_font_3:
					i.font = bold_itatic_24_font_3
	return wb2



############data summary######

###rawdata qc summary
#1.qcsummary(成功)

def rawdata_Q30(inputpath,samplename):
	# make the temp dir
	temp_dir = os.path.join(inputpath, 'tempfile')
	if not os.path.exists(temp_dir):
		os.mkdir(temp_dir)
	tempqc_dir = os.path.join(temp_dir, 'QC')
	if not os.path.exists(tempqc_dir):
		os.mkdir(tempqc_dir)
	qcdir = os.path.join(inputpath, '8Fastqc')
	samplelistdir = os.path.join(qcdir,samplename+'.stat')
	try:
		qcdata = pd.read_table(samplelistdir, nrows=11, sep=' ', names=['reads', 'infor1', 'infor2', 'infor3'])
		qcsummary_re = pd.DataFrame()
		qcsummary_re.loc[0, 'sampleid'] = samplename
		qcsummary_re.loc[0, 'ReadNum'] = qcdata.loc[1, "infor1"].split('\t')[0]
		qcsummary_re.loc[0, 'BaseNum'] = qcdata.loc[1, "infor2"].split('\t')[0]
		qcsummary_re.loc[0, 'Q10'] = qcdata.loc[8, "infor3"]
		qcsummary_re.loc[0, 'Q20'] = qcdata.loc[9, "infor3"]
		qcsummary_re.loc[0, 'Q30'] = qcdata.loc[10, "infor3"]
		qcsummary_re.loc[0, 'GC'] = qcdata.loc[2, "infor1"].split('\t')[0]
		outputname = os.path.join(tempqc_dir,samplename+'_rawQ30.txt')
		qcsummary_re.to_csv(outputname, index=False, header=True, encoding='gbk', sep='\t')
	except:
		print(sampleid + ' is null!')



#2.获得mapping rate summary针对flagstat的结果(成功)
def mapping_rate(inputpath,samplename):
	# make the temp dir
	temp_dir = os.path.join(inputpath, 'tempfile')
	if not os.path.exists(temp_dir):
		os.mkdir(temp_dir)
	tempqc_dir = os.path.join(temp_dir, 'QC')
	if not os.path.exists(tempqc_dir):
		os.mkdir(tempqc_dir)
	ontargetdir = os.path.join(inputpath, '9Ontarget')
	samplelistdir=os.path.join(ontargetdir,samplename+'_clean.abra2.flagstat')
	try:
		samplemapping = pd.read_table(samplelistdir, sep='+', names=['reads', 'infor1', 'infor2'])
		ontarget_re = pd.DataFrame()
		ontarget_re.loc[0, 'sampleid'] = samplename
		ontarget_re.loc[0, 'Read Num'] = samplemapping.loc[0, 'reads']
		ontarget_re.loc[0, 'mapped'] = samplemapping.loc[4, 'reads']
		df = samplemapping["infor1"].str.split('(', expand=True)
		mappingrate = df.loc[4, 1].split(' ')[0]
		ontarget_re.loc[0, 'Mapping Ratio (%)'] = mappingrate
		outputname = os.path.join(tempqc_dir, samplename + '_mappingrate.txt')
		ontarget_re.to_csv(outputname, index=False, header=True, encoding='gbk', sep='\t')
	except:
		print(samplename + ' is null!')


#3.对ontarget rate summary,
#perl /cgdata/guocq/script/wf_idtwes/qualimap_extract.pl -d ${inputpath} bamqc.txt
#ontargetcmd='perl  /cgdata/guocq/script/wf_idtwes/qualimap_extract.pl '+' -d '+inputpath+' '+summarydir+'/'+'table3_ontarget.txt'
#os.system(ontargetcmd)
#inputpath=/cgdata/liuxiangqiong/work62pancancer/pipelinetest-lane3
#perl qualimap_extract.pl -d  ${inputpath} -t p table3_ontarget_rate_1.txt

def ontarget_insertsize(inputpath,samplename):
	laneid = inputpath.split('/')[-1].split('-')[-1]
	# make the temp dir
	temp_dir = os.path.join(inputpath, 'tempfile')
	if not os.path.exists(temp_dir):
		os.mkdir(temp_dir)
	tempqc_dir = os.path.join(temp_dir, 'QC')
	if not os.path.exists(tempqc_dir):
		os.mkdir(tempqc_dir)
	bamdir=os.path.join('/cgdata/pancancer/analyse/',laneid+'/bamfile')
	outputname=os.path.join(tempqc_dir,samplename+'_ontarget1.txt')
	ontargetcmd='perl  /cgdata/liuxiangqiong/work62pancancer/pipeline/v0/prim_umi_v0/qualimap_extract.pl '+' -d '+bamdir+' '+'-t p '+outputname
	os.system(ontargetcmd)
	ontarget = pd.read_table(outputname, sep='\t', header=0)
	ontarget['sampleid'] = ontarget['Sample'].str.split("_", expand=True)[0]
	del ontarget['Sample']
	ontarget_sample=ontarget[ontarget['sampleid']==samplename]
	outputname1 = os.path.join(tempqc_dir, samplename + '_ontarget.txt')
	ontarget_sample.to_csv(outputname1, index=False, header=True, encoding='gbk', sep='\t')
	os.remove(outputname)



###4.获得每个探针的测序深度
def coverage_uniform_samtools(inputpath,samplename):
	# make the temp dir
	temp_dir = os.path.join(inputpath, 'tempfile')
	if not os.path.exists(temp_dir):
		os.mkdir(temp_dir)
	tempqc_dir = os.path.join(temp_dir, 'QC')
	if not os.path.exists(tempqc_dir):
		os.mkdir(tempqc_dir)
	#读入靶区间
	targetdir = '/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/NanOnco_Plus_Panel_v2.0_Covered_b37_cg.parY2X.sort.bed'
	targetdata = pd.read_table(targetdir, sep='\t', header=None, names=['chr', 'start', 'end', 'gene', 'infor1', 'strand'])
	#读入每个靶标区间的测序深度
	coveragedir = os.path.join(inputpath, '10Coverage')
	covdir = os.path.join(coveragedir, samplename + '.cov.samtools_coverage.txt')
	if os.path.getsize(covdir) != 1:
		coverdata_tumor = pd.read_table(covdir, header=0, sep='\t', low_memory=False)
		coverdata_tumor['start'] = coverdata_tumor['start'].astype('int')
		coverdata_tumor['end'] = coverdata_tumor['end'].astype('int')
		tumor_covsum = pd.merge(targetdata, coverdata_tumor, on=['chr', 'start', 'end', 'gene', 'strand'],how='left')
		outputname1 = os.path.join(tempqc_dir,samplename+'_target_coverage.txt')
		tumor_covsum.to_csv(outputname1, index=False, header=True, encoding='gbk', sep='\t')
	# 计算每个样本的平均有效测序深度
	# 读入每个碱基位点测序深度和均一性
	coveragedir = os.path.join(inputpath, '10Coverage')
	unfilterreadsdir = os.path.join(coveragedir, samplename + '.cov.samtools.txt')
	if (os.path.exists(unfilterreadsdir)) and (os.path.getsize(unfilterreadsdir) != 0):
		inputdata1 = pd.read_table(unfilterreadsdir, header=None, sep='\t', low_memory=False)
		inputdata1.columns = ['chrom', 'start', 'readcounts']
		##获得样本的平均测序深度
		# T_mean = ontarget[ontarget['sampleid'] == sampleid].reset_index(drop=True).loc[0, 'T_Mean']
		T_mean = int(inputdata1['readcounts'].sum() / len(inputdata1))
		# 计算覆盖均一性（大于平均深度20%的碱基位点占目标区域碱基位点总数的比例）
		Coverge_uniform = len(inputdata1[inputdata1['readcounts'] > 0.2 * T_mean]) / len(inputdata1)
		averagedata = pd.DataFrame()
		averagedata.loc[0, 'sampleid'] = samplename
		averagedata.loc[0, 'Unique_average_depth'] = round(T_mean,4)
		averagedata.loc[0, 'coverge_uniform']=round(Coverge_uniform,4)
		outputname2 = os.path.join(tempqc_dir, samplename + '_avarageDepth_coverage_uniform.txt')
		averagedata.to_csv(outputname2, index=False, header=True, encoding='gbk', sep='\t')



####输出原始的测序深度
#本panel大小2.6M，编码区大小1.39M
#总体质控合格的判断
#1)总的平均原始测序深度血液>=5000X;组织>=1000X;白细胞对照>=200X
#2)有效测序深度，血液>=1000X;组织>=500X
#2)血液insert size （median）<=200;组织insert size(median)>=100
#3)覆盖均一性>=80%
#4)Q30>=80%
####补充计算覆盖均一性
#大于平均深度20%的碱基位点占目标区域碱基位点总数的比例

####修改采用samtools计算平均测序深度时候

def QCreprot2_samtools(inputpath,samplename):
	laneid = inputpath.split('/')[-1].split('-')[-1]
	temp_dir = os.path.join(inputpath, 'tempfile')
	tempqc_dir = os.path.join(temp_dir, 'QC')
	#inputpath = '/cgdata/liuxiangqiong/work62pancancer/pipelinetest-lane3'
	# 本panel大小2.6M
	panelsize = 2600000
	# 读入基础质控
	rawdatadir=os.path.join(tempqc_dir,samplename+'_rawQ30.txt')
	QCdata = pd.read_table(rawdatadir, sep='\t', header=0)
	###读入insert size
	insertsizedir = os.path.join(tempqc_dir, samplename + '_ontarget.txt')
	ontarget = pd.read_table(insertsizedir, sep='\t', header=0)
	# 读入平均有效测序深度和覆盖均一性
	coveragedir = os.path.join(tempqc_dir,  samplename + '_avarageDepth_coverage_uniform.txt')
	coveragedata = pd.read_table(coveragedir, header=0, sep='\t')
	##获得样本的有效平均测序深度
	T_mean = coveragedata.loc[0, 'Unique_average_depth']
	# 获得均一性
	Coverge_uniform = coveragedata.loc[0, 'coverge_uniform']
	# 获得insert size(median)
	Insert_size_median = int(ontarget.loc[0, 'Insert_Size'].split(' / ')[1])
	##获得Q30
	Q30 = QCdata.loc[0, 'Q30']
	##获得原始数据量
	rawdata_base = 2 * 150 * (QCdata.loc[0, 'ReadNum'])
	# 获得原始测序深度
	raw_depth = round(rawdata_base / panelsize, 4)
	# 总体质控合格的判断
	QC_report=pd.DataFrame()
	sampletype = samplename[-2:]
	if sampletype == 'TT':
		if (raw_depth >= 1000) & (Insert_size_median >= 100) & (Coverge_uniform >= 0.8) & (
				(float(Q30.split('%')[0]) / 100) >= 0.8) & (T_mean>=500):
			label = 'Success'
		else:
			label = 'Fail'
	elif sampletype == 'CT':
		if (raw_depth >= 5000) & (Insert_size_median <= 200) & (Coverge_uniform >= 0.8) & (
				(float(Q30.split('%')[0]) / 100) >= 0.8)& (T_mean>=1000):
			label = 'Success'
		else:
			label = 'Fail'
	elif sampletype == 'CN':
		if raw_depth >= 200:
			label = 'Success'
		else:
			label = 'Fail'
	QC_report.loc[0, 'sampleid'] = samplename
	QC_report.loc[0, 'Total_base'] = rawdata_base
	QC_report.loc[0, 'Total_average_depth'] = round(raw_depth)
	QC_report.loc[0, 'Unique_average_depth'] = T_mean
	QC_report.loc[0, 'insert_size'] = Insert_size_median
	QC_report.loc[0, 'coverge_uniform'] = str(round(Coverge_uniform * 100, 2)) + '%'
	QC_report.loc[0, 'Q30'] = Q30
	QC_report.loc[0, 'QC_overall'] = label
	print(QC_report)
	outputname = os.path.join(tempqc_dir, samplename + '_qc_report.txt')
	QC_report.to_csv(outputname, index=False, header=True, encoding='gbk', sep='\t')
	# 输出疾病样本的qc报告。那么需要改变一下样本名，然后再输出
	if (sampletype == 'TT') or (sampletype == 'CT'):
		samplenewid = samplename[:-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, samplenewid)
		if not os.path.exists(sample_dir):
			os.mkdir(sample_dir)
		QC_report.loc[0, 'sampleid'] = samplenewid
		QC_report1 = QC_report[
			['sampleid', 'Total_average_depth', 'insert_size', 'coverge_uniform', 'Q30', 'QC_overall']]
		outputfile1 = os.path.join(sample_dir, laneid + '-' + samplenewid + '.qc1.xlsx')
		writer = pd.ExcelWriter(outputfile1)
		QC_report1.to_excel(writer, sheet_name='QC', index=False)
		writer.save()
		writer.close()
		if sampletype == 'TT':
			sampleQCdata = excel_font_color(outputfile1, "B", 2, 1000, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "C", 2, 100, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "D", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "E", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "F", 3, "Fail", "red")
		elif sampletype == 'CT':
			sampleQCdata = excel_font_color(outputfile1, "B", 2, 5000, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "C", 1, 200, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "D", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "E", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "F", 3, "Fail", "red")
		elif sampletype == 'CN':
			sampleQCdata = excel_font_color(outputfile1, "B", 2, 200, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "C", 1, 180, "red")
			sampleQCdata = excel_font_color(sampleQCdata, "D", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "E", 2, 0.8, "red", percentage=True)
			sampleQCdata = excel_font_color(sampleQCdata, "F", 3, "Fail", "red")
		outputfile2 = os.path.join(sample_dir, laneid + '-' + samplenewid + '.qc.xlsx')
		sampleQCdata.save(outputfile2)
		os.remove(outputfile1)







#主要的执行命令
def qcrun(inputpath,samplename):
	print('step1.calculate the Q30')
	rawdata_Q30(inputpath,samplename)
	print('step2.mapping rate')
	mapping_rate(inputpath,samplename)
	print('step3.insertsize')
	ontarget_insertsize(inputpath,samplename)
	print('step4.coverage uniform')
	coverage_uniform_samtools(inputpath,samplename)
	print('step5.QCreport')
	QCreprot2_samtools(inputpath,samplename)



if __name__=='__main__':
	parser = argparse.ArgumentParser(description='for the QC')
	parser.add_argument('-i', '--inputpath', type=str, help='the path of lane')
	parser.add_argument('-s', '--samplename', type=str, help='all of the sample')
	args = parser.parse_args()
	Inputpath = args.inputpath
	Samplename = args.samplename
	qcrun(Inputpath,Samplename)