外显子（wes）panel数据分析OncoKB注释

作者，Evil Genius

接下来我们进行第二步，对外显子突变信息进行OncoKB注释，为第三步用药指导做充分的准备。我们进行的分析从annovar注释后的结果开始，这里注意很多的细节就可以了。

这里的脚本是python写的全自动化脚本，首先我们需要定义需要的处理函数，这其中包括的函数如下：

函数mutation_vf：计算突变的频率
函数df_aachange : vcf_filter文件的aachange列按p.进行去重，保留NM号较小的值
函数df_geneDetail:vcf_filter文件的geneDetail列按c.进行去重，保留NM号较小的值
函数get_oncokb_maf_file: annovar注释后的txt文件转换成OncoKB注释的矩阵文件。

完整的定义过程及自动化脚本放在了最后,我们先来解析脚本

import os
import re
# import time
from openpyxl import load_workbook
id = sample_name
input_file = os.path.join(current_path, '检测结果', id+'.filter.vcf.hg19_multianno.txt')
output_file = os.path.join(current_path, '检测结果', id+'.oncokb.maf.txt')
get_oncokb_maf_file(input_file,output_file)
oline = '\t'.join([id,'NSCLC'])
clinical_info_file.write(oline+'\n')

这一步中输入annovar的注释结果（txt格式），示例如下：

输出我们需要的OncoKB注释需要的矩阵文件，如下图

OncoKB在先注释，这里展示多脚本的注释方法，其实就是循环，每个id为每个样本。

for id in ids:
    # time.sleep(5)
    print('python3 MafAnnotator.py -i 检测结果/' + \
        id + '.oncokb.maf.txt -o 检测结果/' + id + '.oncokb.txt -c clinical_info.txt -b "84a75c6f-7daf-4fee-86ce-ae3f12f31413" -q hgvsg')

    try:
        os.system('python3 MafAnnotator.py -i 检测结果/' + \
            id + '.oncokb.maf.txt -o 检测结果/' + id + '.oncokb.txt -c clinical_info.txt -b 84a75c6f-7daf-4fee-86ce-ae3f12f31413 -q hgvsg')
    except:
        print('''
            
            
            %s样本报错，请注意！


        ''' % id)
        continue

# 处理在线oncokb注释后的文件
    oncokb_file = open(os.path.join(current_path,'检测结果',id+'.oncokb.txt'),'r')
    oncokb_file_tmp = open(os.path.join(current_path,'检测结果',id+'.oncokb.tmp'),'w')

    for line in oncokb_file:
        array = line.split('\t')
        if array[0] == 'NCBI_Build':
            array[11:13] = ['sampleNum','sampleid']
            array_t = array[11:16]
            array_t.reverse()
            array[11:16] = array_t
            oarray = array[:16] + ['sampleRatio','avsnp150','cosmic95','fathmm-MKL_coding_score','fathmm-MKL_coding_pred'] + array[16:]
        else:
            array_t = array[11:16]
            array_t.reverse()
            array[11:16] = array_t
            out_array = array[7:12]
            temp_k = '~'.join(out_array)
            temp_k = 'chr' + temp_k
            # print(temp_k)
            if temp_k in data:
                samples = ';'.join(data[temp_k])+';'
                temp_num = len(data[temp_k])
                array[14:16] = [samples,str(temp_num)+'/'+str(idsSum)]
            else:
                temp_num = 0
                array[14:16] = ['.','0/'+str(idsSum)]
        
            oarray = array[:16] + [str("{:.2f}%".format(temp_num/idsSum*100))] + vcf_info[temp_k] + array[16:]

        oline = '\t'.join(oarray)
        oncokb_file_tmp.write(oline)

    oncokb_file.close()
    oncokb_file_tmp.close()


    os.remove(os.path.join(current_path,'检测结果',id+'.oncokb.txt'))
    os.rename(os.path.join(current_path,'检测结果',id+'.oncokb.tmp'),os.path.join(current_path,'检测结果',id+'.oncokb.txt'))

这样就得到了OncoKB的注释结果，部分截图如下：

完成了这个，我们就完成了clinvar、intervar、cosmic、oncoKB的位点注释，完全足够知到下游的用药指导了。

封装好以后直接全自动出结果即可，下面是示例

#####进入annovar注释的目录下
python3  主脚本即可

完整主脚本如下：

# -*- coding: utf-8 -*-

import os
import re
# import time
from openpyxl import load_workbook

def mutation_vf(info):
    # 计算突变的频率
    arr_t = info.split(':')[1].split(',')
    arr_t = [int(i) for i in arr_t]
    dp = int(info.split(':')[3])
    alt_dp = max(arr_t[1:])
    vf = alt_dp/dp

    return vf

def df_aachange(aachange):
    '''
    description: vcf_filter文件的aachange列按p.进行去重，保留NM号较小的值；
        # aachange = 
        # 'FGFR2:NM_001144916:exon2:c.T212C:p.M71T,FGFR2:NM_001144918:exon3:c.T212C:p.M71T,FGFR2:NM_023029:exon3:c.T290C:p.M97T,
        # FGFR2:NM_001144913:exon4:c.T557C:p.M186T,FGFR2:NM_001144914:exon4:c.T557C:p.M186T,FGFR2:NM_001144915:exon4:c.T290C:p.M97T,
        # FGFR2:NM_001144919:exon4:c.T290C:p.M97T,FGFR2:NM_000141:exon5:c.T557C:p.M186T,FGFR2:NM_001144917:exon5:c.T557C:p.M186T,
        # FGFR2:NM_001320658:exon5:c.T557C:p.M186T,FGFR2:NM_022970:exon5:c.T557C:p.M186T'
    param {aachange：aachange列信息；}
    return {aachange_uni