来源https://github.com/sigven/pcgr
使用注意事项
1.若使用软连接,在该文本中无法识别,所以在mount上docker的镜像时会失去软连接的能力,从而导致数据库无法被mount在镜像中,从而导致无法使用。
代码digging
- pcgr.py为主要接口,主要在接收参数并且构造docker run的命令从而恰当的mount完所有的input/output。
- 接收参数并检查输入输出的文件夹是否存在。
- 构建Logger
- 创建子进程,handle子进程的错误
- 从参数中判断并构造恰当的docker run的命令
- 其中pcgr.py 的 执行流程涉及的py文件。
- pcgr_check_input.py 使用vep对其进行注释
- pcgr_vcfanno.py 用多个数据库对其进行关联注释
- pcgr_summarise.py 总结注释后的结果
- pcgr.R 生成总结性的结果,接收一个通过vcfanno注释后的vcf文件,其中包括了Tier位点,引用文献。
- 从执行的命令中可以看出执行顺序以及产生的文档。
- pcgr_check_input.py
- VEP使用较为复杂参数的注释,其中有一个output的该vcf文件注释过程中所产生的简单html类型的数据统计
- 用sed删除一部分vep进行注释后再vcf中所产生的会干扰后续分析流程的参数
- 使用pcgr_vcfanno.py进行整合注释
pcgr_check_input.py 具体内容
- 接收简单的参数,如pcgr的目录,input_vcf / input_cna
- 验证cna / vcf的完整程度
- 验证是否是合格的pcgr的输入文件,若不是,对其进行校正。
- 用EBIvariation/vcf-validator 对其进行校验是否为正确的vcf
- 检查INFO tags以了解是否会重复注释
- 若该vcf中含有multipe alternative alleles, 用vt工具进行decomposition
- 检查cna文件的列是否合格,以及数据格式
- 检查是否排序并做好索引,否则用bgzip与tabix进行处理
pcgr_vcfanno.py 具体内容
- 接收参数,并使用vcfanno同时对多个vcf文件进行平行注释,主要接收的是查询的vcf,输出的vcf,以及pcgr的根目录,还有所有数据库的指定(作者已经将其统一为了vcf文件)
- 构建一个head文件临时存储head信息,最后再一起整合进入新的vcf文件。
- appedn_to_conf_file专门给输出的配置文件指明用了哪些数据库进行注释,注释的过程如何,属于临时文件,阅后即焚。
- 辅助函数,get_vcf_info_tags与print_vcf_header
- 主要进行注释的脚本
若需要进行数据库增减修改,主要hack入vcfanno这个软件的参数,按要求提供所需要的注释数据库。并修改该文件的入口。后续总结脚本中应该还有接口需要注意。
- 主函数的流程:
pcgr_summarise.py
- 接收参数: 已经进行注释好的vcf文件,以及pcgr的根目录
- 主函数: extend_vcf_annotations
- 辅助函数:
- 将3字母氨基酸名称转化为单字母
- 通过dbnsfp对突变效果进行预测并注释
- 翻译蛋白质信息,转化为特定的columns,加入表头
- 给定transcript ID ,转化其为多行基因相关信息。
-
主函数的流程:
- 加载数据库
涉及数据库,若修改,需要注意。
2. 在input vcf的当前目录下创建注释后的vcf
3. 读入记录vep和pcgr的相应表头信息的表格。
4. 从vcf的头部信息中分析vep注释后的头部信息,并且初始化两个index与该数据库名称的(相互关系)字典。
5. 利用dbnsfp对突变效果进行预测(主要是利用了SIFT、PolyPhen2等评价数字)
6. 加入EFFECT_PREDICITONS的INFO行
7. 将3中读入的INFO头部信息加入vcf INFO。
8. 逐行处理variant,若改行未被VER注释则跳过,通过辅助函数,从vep的注释INFO头部扩展出多个蛋白、基因相关的INFO头部。
9. 设置该INFO头部下的值。共计5批tags INFO
- VEP的信息,即Consequence annotations from Ensembl VEP
- dbNSFP的注释信息
- 原始input vcf本就存在的info
- pcgr中癌症相关的基因注释
- pcgr中蛋白相关的基因注释,如(hotspots, features etc.)
10. 处理完同一染色体的突变后,统一写入output vcf,以染色体为单位写入一批突变位点注释信息。
11. 压缩并建立索引
pcgrr2 R语言脚本的内部结构
共计17个函数
按照调用顺序阐述
1. generate_pcg_report
- 初始化列表和输出的中间文件
- 调用
get_calls
读取vcf文件,储存入对象sample_calls- 使用
generate_report_data
转换sample_calls为report_data- 检查report_data并输出其中的部分matrix到中间文件
- 若需要输出成report,则将report的部分信息传输给R的库
rmarkdown
的render函数,渲染出一个完整的html页面。
2. get_calls
- 使用BioConductor的包
BSgenome.Hsapiens.UCSC.hg19
来进行读取vcf,转换vcf- 其中使用了
VariantAnnotation::called
进行了一步筛选,尤其注意的是,其中的过滤为默认选项,并且测试中,将12000+的vcf过滤掉后只剩下600+的位点,其中丢弃了较为重要的具有生物学意义的位点。若将该行注释掉并不影响后续代码。- 若无数据行,则需要把应有的列加上去。
- 使用
postprocess_vranges_info
对其中的区间信息进行预处理- 加入pfam的数据库关联
add_pfam_domain_links
- 加入swissport数据库关联
add_swissprot_feature_descriptions
- 通过
dplyr:: left_join
将现有的vcf与其VAR_ID列与作者提供的pcgr_data进行关联扩展。- 同上的对ENTREZ_ID与gene进行关联扩展。
- 进行一些命名上的修改,包括去除chr,重命名
- 关联扩展CLINVAR,其中包括去除了原数据中的特殊字符。
- 关联扩展了各个数据的links ,通过
annotate_variant_link
- return big data frame
3.generate_report_data
- 接收一个R中的data.frame,初始化所有的参数。
- 挑选出coding var、noncoding var、SNV、INDEL,各自单独储存
- 选出除去MT染色体上的所有SNV位点作为signature_call_set并用
signature_contributions_single_sample
对其进行计算权重???并且作为一个衡量标准,由于其中其中每一个signature代表一类疾病,通过pcgr_data$signatures_aetiologies的Signature_ID进行关联,判断得知更有可能为哪一类疾病。结束后,将其余原始的权重的矩阵合并成一个list,储存下来。- Tier 1: 使用
get_clinical_associations_civic_cbmdb
对coding var得到举证进行对civic的关联,将输出储存如clinical_evidence_items_tiers1A/1B/1C矩阵中,(分类标准见get_clinical_associations_civic_cbmdb
)。合并后作为一个完整的variants_tier1矩阵- Tier 2: 在coding var的基础上,将
INTOGEN_DRIVER_MUT & CANCER_MUTATION_HOTSPOT & OTHER_DISEASE_DOCM
为空的位点去除后,再去除与Tier1重复的位点,后剩下的位点。- Tier3: 在coding var得到基础上,选择CANCER_CENSUS_SOMATIC非空,或者ONCOGENE为真,TUMOR_SUPPRESSOR为真的位点,去重后,作为Tier3
- Tier4: 去除前3个tier后的剩下的所有的coding位点
- Tier5: 所有noncoding var
- 将5个tier所得到的位点通过
tier_to_maf
整合成maf,用generate_biomarker_tsv
将tier1部分输出为biomarker_tsv。- 将用到的大部分东西整合成一个巨大的list,report_data,并返回。
4. postprocess_vranges_info
- 从一个特殊的对象中提取出chr,pos,start,alt,用_将其串联在一起作为新的一列,即VAR_ID
5. add_pfam_domain_links
- 选取出DOMAINS列不为空的行,存入新的对象
- 去重,去除冗余信息,将pcgr_data中的关联信息用left_join加入该对象。
- 再次left_join将新对象和旧对象合在一起。
6.add_swissprot_feature_descriptions
- 同上
7. annotate_variant_link
- 将各个数据库的来源地址附加上到矩阵中,不同的数据库最后添加上去的列名都不一样,用了5次判断语句进行选择。
8. signature_contributions_single_sample
- 主要为
deconstructSigs
函数调用deconstructSigs::mut.to.sigs.input
将输入的数据提取出碱基替换情况,生成行为sample,列为每3个碱基中的碱基替换情况,为cosmic中分的mutation signature(96=644),即嘧啶的替换(6)左边方向(4)右边方向的(4)deconstructSigs::whichSignatures
从中结合自身包中含有的signatures.cosmic
的特征权重值,将其与现有的碱基替换情况进行计算,以确定能用来重建突变profile最适合的signatures.cosmic
中提供的30个signature的权重值。- 提取出权重的表,并返回。
9. get_clinical_associations_civic_cbmdb
- 预处理pcgr_data
- 其中包括三种mapping方式
- exact方式,先挑出CIVIC_ID列不为空的dataframe。
- 若CIVIC_ID中含有EID,通过CIVIC_ID与
pcgr_data$civic_biomarkers
中的evidence_id进行关联join,否则与pcgr_data$civic_biomarkers
中的civic_id进行关联join。- vcf_data_df_cdmdb为vcf_data_df中CIVIC_ID和CBMDB_ID均不为空的进行初始化
- 通过CIVIC_ID与
pcgr_data$cbmdb_biomarkers
中的CBMDB_ID进行关联join
- approximate方式,其中又可以细分为codon or exon-level,先挑出CIVIC_ID_2列不为空的dataframe。
- EID部分的关联类似。后把
pcgr_data$civic_biomarkers
中mapping_category部分等于gene的部分筛掉。- 后续步骤主要为过滤掉非空的CODON或者EXON列,并且进行一定的修改以符合后续的分析。
- 最后从pcgr_data中并加入新的列名,并且其作为切片标准的一部分。从merge后的大表中选取特定列作为输出并返回。
10.tier_to_maf
- 主要作用时给特定的CONSEQUENCE下的位点进行一个位点的分类,分类完以后还有一个固定的字符串顺序以进行给染色体的排序。
修改想法与记录。
17.06.20
- 注释掉get_calls中筛选语句,从而保留下大部分的位点,尤其是T790M等重要位点。