pcgr使用记录(主要为源码分析)

2017-06-13  本文已影响0人  栽生物坑里的信息汪

来源https://github.com/sigven/pcgr

使用注意事项

1.若使用软连接，在该文本中无法识别，所以在mount上docker的镜像时会失去软连接的能力，从而导致数据库无法被mount在镜像中，从而导致无法使用。

代码digging

  1. pcgr.py为主要接口，主要在接收参数并且构造docker run的命令从而恰当的mount完所有的input/output。
  1. 其中pcgr.py 的 执行流程涉及的py文件。
  2. pcgr_check_input.py 使用vep对其进行注释
  3. pcgr_vcfanno.py 用多个数据库对其进行关联注释
  4. pcgr_summarise.py 总结注释后的结果
  5. pcgr.R 生成总结性的结果，接收一个通过vcfanno注释后的vcf文件，其中包括了Tier位点，引用文献。
  6. 从执行的命令中可以看出执行顺序以及产生的文档。
  1. pcgr_check_input.py
  1. VEP使用较为复杂参数的注释，其中有一个output的该vcf文件注释过程中所产生的简单html类型的数据统计
  2. 用sed删除一部分vep进行注释后再vcf中所产生的会干扰后续分析流程的参数
  3. 使用pcgr_vcfanno.py进行整合注释

pcgr_check_input.py 具体内容

  1. 接收简单的参数，如pcgr的目录，input_vcf / input_cna
  2. 验证cna / vcf的完整程度
  3. 验证是否是合格的pcgr的输入文件，若不是，对其进行校正。
    • 用EBIvariation/vcf-validator 对其进行校验是否为正确的vcf
    • 检查INFO tags以了解是否会重复注释
    • 若该vcf中含有multipe alternative alleles， 用vt工具进行decomposition
    • 检查cna文件的列是否合格，以及数据格式
    • 检查是否排序并做好索引，否则用bgzip与tabix进行处理

pcgr_vcfanno.py 具体内容

  1. 接收参数，并使用vcfanno同时对多个vcf文件进行平行注释，主要接收的是查询的vcf，输出的vcf，以及pcgr的根目录，还有所有数据库的指定（作者已经将其统一为了vcf文件）
  2. 构建一个head文件临时存储head信息，最后再一起整合进入新的vcf文件。
  3. appedn_to_conf_file专门给输出的配置文件指明用了哪些数据库进行注释，注释的过程如何，属于临时文件，阅后即焚。
  4. 辅助函数，get_vcf_info_tags与print_vcf_header
  5. 主要进行注释的脚本

若需要进行数据库增减修改，主要hack入vcfanno这个软件的参数，按要求提供所需要的注释数据库。并修改该文件的入口。后续总结脚本中应该还有接口需要注意。

  1. 主函数的流程：

pcgr_summarise.py

  1. 接收参数： 已经进行注释好的vcf文件，以及pcgr的根目录
  2. 主函数： extend_vcf_annotations
  3. 辅助函数：

  1. 主函数的流程：

    1. 加载数据库

涉及数据库，若修改，需要注意。

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
  1. VEP的信息，即Consequence annotations from Ensembl VEP
  2. dbNSFP的注释信息
  3. 原始input vcf本就存在的info
  4. pcgr中癌症相关的基因注释
  5. 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

  1. 注释掉get_calls中筛选语句，从而保留下大部分的位点，尤其是T790M等重要位点。
上一篇 下一篇

猜你喜欢

热点阅读