实验记录 | 梳理代码2

凡事预则立，不预则废。我就不相信我多下一点功夫会搞不定它。我一定要把它弄清楚。
刚刚看到github上关于这个流程的讨论，对于作者佩服的不行。
觉得只有更加的努力，才能在这个领域有自己的一席之地。想想看已经有学弟比自己厉害了。既然学弟可以顺利的跑下来，那么我想，我也可以。我想成为这个领域最厉害的那20%。进阶之路真的是永无止境啊。
这篇文章，我打算梳理一下github上，ORBC pipeline的框架。都主要有哪一些的文件。并且每一个文件的意思是什么？
那个学弟之所以可以与作者“对话”，主要原因还是细致的看了源代码。既然源代码在那里，大家的学习资料都是一样的，而且我也有这么一些编程基础。为什么他可以，我不可以呢？

参考链接：https://github.com/tianshilu/QBRC-Somatic-Pipeline/issues/7

由于我现在电脑中，还在从百度网盘上下载文件，已经下载到71%，是实在不想前功尽弃。于是，还是在window平台上进行梳理。

E:.
├─disambiguate_pipeline
│ └─conda_env
│ ├─bin
│ ├─conda-meta
│ ├─include
│ │ └─bamtools
│ │ ├─api
│ │ │ └─algorithms
│ │ └─shared
│ ├─lib
│ │ └─pkgconfig
│ ├─share
│ │ ├─info
│ │ └─licenses
│ │ ├─gcc-libs
│ │ └─libstdc++
│ └─x86_64-conda_cos6-linux-gnu
│ └─sysroot
│ └─lib
├─example
│ ├─cnv
│ ├─cnv_summarization
│ ├─example_dataset
│ │ ├─example_output
│ │ └─sequencing
│ └─filter
│ └─each
└─somatic_script
一级目录下，还有以下独立的文件

• perl
  • cnv.pl
  • evolution.pl （未说明）
  • job_cnv.pl
  • job_somatic.pl
  • somatic.pl
• R
  • design_summarize_cnv.R （未说明）
  • filter.R
  • plot_cnv.R （未说明）
  • summarize_cnv.R

还有两版的readme文件，是对流程的一个说明与补充。
我们一点点的开始分析。

1。README

一共有两个README文档。其中README.md之前已经分析过。主要的内容是介绍每一个文件的使用规则。
而从另一个README文档中得到的信息主要是：

• 该分析流程适用的物种：人类和小鼠
• matched normal sample
  有匹配的正常样本的mutation calling比较严格。而没有匹配的正常样本的mutation calling则相对宽松。
• 其中包含有一个evolution.pl的代码文件。这个文件根据来自于一个患者的一系列的样本的点突变以及拷贝变异的结果，来推测进化关系。这个文件是后来的SClineagar的雏形吗？值得一提。
• 也包含着更新的日志，在何时添加了哪一个变量，对于系统整体的功能丰富在什么地方。

这是我对于这个文档的理解。

2。4个perl文件。

并非是浅尝辄止的，而是要把文件打开来看看里面是什么？怎样的调用程序的，他的结构是什么样的？参数，设定的阈值？等等信息。只有把这个层次弄清楚了，才能很好的驾驭好这个流程。

（1）somatic.pl

这张图片，展示的就是somatic.pl的处理的整个过程。也是我们这次pipeline最核心的一个部分。
其输入文件是分别来源于正常组织和肿瘤组织的fastq.gz文档。输出文档就是点突变的mutation calling的结果。
我们先打开看一下测试数据集的mutation calling的结果是什么样子的（文件存放在\example\example_dataset\example_output这个文件路径下了）。
一共有两个文档，分别命名为：（1）germline_mutations_hg38.txt（2）somatic_mutations_hg38.txt
首先需要弄明白的是，germline突变与somatic突变的区别。
按照我的理解，简单的总结一下：
（1）germline和somatic突变都属于点突变的范畴。但是产生的原因不一样，造成的结果也不一样。
（2）germline是生殖细胞的突变，主要是在胚胎发育的过程中，生殖细胞（精子，卵细胞）的基因突变所致。这个突变所导致的结果即是，该个体所有细胞都包含有这种突变。而只有通过与其他正常的个体参照，才能发现这类突变。并且其VAF是很大的。
（3）somatic则是体细胞随机突变的结果。可以在任意体细胞中发生，但是相对的，其VAF非常小（因为是随机的过程，且范围在全身）。

弄清楚这一点之后，我们来看一下个自己的calling结果文档是什么样子的。
（1）germline_mutations_hg38.txt

Chr Start End Ref Alt Caller Normal_ref Normal_alt Tumor_ref Tumor_alt Func.refGene Gene.refGene ExonicFunc.refGene AAChange.refGene SIFT_pred Polyphen2_HVAR_pred cosmic70 esp6500siv2_all ExAC_ALL X1000g2015aug_all

==========================================================================
chrM 3394 3394 T C lofreq 528 4 528 4 downstream RNR2 . . . . 0 0 0 0
chrM 3460 3460 G A lofreq 394 5 394 5 downstream RNR2 . . . . 0 0 0 0
chrM 3474 3474 C A lofreq 443 4 443 4 downstream RNR2 . . . . 0 0 0 0
chrM 3670 3670 G T lofreq 436 3 436 3 downstream RNR2 . . . . 0 0 0 0
chrM 3892 3892 A AC lofreq 539 5 539 5 downstream RNR2 . . . . 0 0 0 0
chrM 4048 4048 G A lofreq 1 82 1 82 downstream RNR2 . . . . 0 0 0 0
chrM 4071 4071 C T lofreq 0 260 0 260 downstream RNR2 . . . . 0 0 0 0
chrM 4164 4164 A G lofreq 2 277 2 277 downstream RNR2 . . . . 0 0 0 0
chrM 4769 4769 A G lofreq 2 148 2 148 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 5351 5351 A G lofreq 1 128 1 128 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 5460 5460 G A lofreq 1 222 1 222 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 6455 6455 C T lofreq 1 459 1 459 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 6680 6680 T C lofreq 0 516 0 516 intergenic RNR2;NONE . . . . 0 0 0 0
……

（2）somatic_mutations_hg38.txt

Chr Start End Ref Alt Caller Normal_ref Normal_alt Tumor_ref Tumor_alt Func.refGene Gene.refGene ExonicFunc.refGene AAChange.refGene SIFT_pred Polyphen2_HVAR_pred cosmic70 esp6500siv2_all ExAC_ALL X1000g2015aug_all

=======================================================
chrM 3394 3394 T C lofreq 528 4 528 4 downstream RNR2 . . . . 0 0 0 0
chrM 3460 3460 G A lofreq 394 5 394 5 downstream RNR2 . . . . 0 0 0 0
chrM 3474 3474 C A lofreq 443 4 443 4 downstream RNR2 . . . . 0 0 0 0
chrM 3670 3670 G T lofreq 436 3 436 3 downstream RNR2 . . . . 0 0 0 0
chrM 3892 3892 A AC lofreq 539 5 539 5 downstream RNR2 . . . . 0 0 0 0
chrM 7884 7884 T A lofreq 411 4 411 4 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 9531 9531 A AC lofreq 376 10 376 10 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 9627 9627 G A lofreq 428 221 428 221 intergenic RNR2;NONE . . . . 0 0 0 0
chrM 9778 9778 G T lofreq 561 20 561 20 intergenic RNR2;NONE . . . . 0 0 0 0
……

这两个文件的列名是一致的。那我们就一起来分析。

• Chr：突变所在的染色体
• Start：染色体上的起始位点
• End：染色体上的终止位点
• Ref：参考基因组在该位点上的碱基
• Alt：被替换的碱基
• Caller：突变检测的工具
• Normal_ref/Normal_alt/Tumor_ref/Tumor_alt：这个我没怎么懂，以及下方的数字
• Func.refGene/Gene.refGene/ExonicFunc.refGene/AAChange.refGene/SIFT_pred Polyphen2_HVAR_pred/cosmic70/esp6500siv2_all/ExAC_ALL/X1000g2015aug_all：这一堆，我觉得像是annovar注释信息。

所以从结果文件上，我们可以明白。somatic.pl这个程序就是把输入文件中的突变位点找到。并对其进行注释。只不过，germline/somatic 突变的calling机制会有所区别。
以上，是我现在的理解，只是还有一点存疑：

tumor/normal 到底有啥用呢？——>我或许可以问一下作者。——>发邮件尝试问过，现在没有回复。估计不再继续维护了。

Dear Professors@wtwt5237:
Recently,I have been learning this pipeline,and hope to transplant it in our own data.However,I have some questions about the input data in the somatic.pl.Why you set both the normal and tumor samples at the same time?Does it means that compare the tumor sample with the normal sample,and turn out the mutations in tumor samples against the normal?
I also notice your annoucement that"For tumor-only calling, put “NA NA” in the slots of the normal samples. Results will be written to germline files",Or maybe we can use the tumor only to call the germline files,while use the normal only to call what?
I can not understand the pair of normal and tumor samples,can how to define?The cells that come from the normal and tumor tissue from one patient?or the cells form the normal and camer patients respectively?
In other words,If I want to call the normal person’s somatic mutations in one particular tissues to traces their development lineage,How can I input my files?
Hope for your suggestions!Thank you for your nice work!
Best regrads,
Xiu

我们继续研读这部分的代码。
这个文件是由perl写成的，一共是539行。