肺癌分析插件Can_28修复问题

# 肺癌分析插件Can_28修复问题
## 一、报错首先确定下网页上

（如http://172.24.3.5/report/271 ），barcode对应到的样品名是否拆分出数据，可能是tag用错了或者plan设错了。（一般我们设置好barcode系统后，没有变更过就没问题）

![](/media/202301/2023-01-09_140525&8QPLuEqUwgCjS2yhTIrs.png)

## 二、插件总体思路

对bam文件生成mpileup格式，并分析mpileup结果，找出所有候选变异位点。之后应用1000人的阴性数据库进行变异检测。
### 1、捕获区域为：/results/plugins/Can_28/design/CHP2.20131001.noh.bed，共28个区间，如下：

```
chr1    115256482   115256546   Can_NRAS_a  .   GENE_ID=NRAS
chr1    115258721   115258811   Can_NRAS_b  .   GENE_ID=NRAS
chr1    162748338   162748411   Can_DDR2    .   GENE_ID=DDR2
chr3    178936051   178936111   Can_PIK3CA_a    .   GENE_ID=PIK3CA
chr3    178952026   178952090   Can_PIK3CA_b    .   GENE_ID=PIK3CA
chr7    55241645    55241725    Can_EGFR_18 .   GENE_ID=EGFR
chr7    55242445    55242510    Can_EGFR_19 .   GENE_ID=EGFR
chr7    55248988    55249083    Can_EGFR_20 .   GENE_ID=EGFR
chr7    55259485    55259549    Can_EGFR_21 .   GENE_ID=EGFR
chr7    140453083   140453165   Can_BRAF_15 .   GENE_ID=BRAF
chr7    140481395   140481473   Can_BRAF_11 .   GENE_ID=BRAF
chr10   89717623    89717710    Can_PTEN_a  .   GENE_ID=PTEN
chr12   25380240    25380307    Can_KRAS_3  .   GENE_ID=KRAS
chr12   25398237    25398310    Can_KRAS_2  .   GENE_ID=KRAS
chr14   105246513   105246585   Can_AKT1_4  .   GENE_ID=AKT1
chr17   37880971    37881058    Can_ERBB2_20    .   GENE_ID=ERBB2
chr15   66727443    66727527    Can_MEK1_2  .   GENE_ID=MEK1
chr17   7573984 7574075 Can_TP53_2  .   GENE_ID=TP53
chr17   7576950 7577032 Can_TP53_3  .   GENE_ID=TP53
chr17   7577082 7577173 Can_TP53_4  .   GENE_ID=TP53
chr17   7577511 7577595 Can_TP53_5  .   GENE_ID=TP53
chr17   7578144 7578226 Can_TP53_6  .   GENE_ID=TP53
chr17   7578385 7578474 Can_TP53_7  .   GENE_ID=TP53
chr17   7579290 7579381 Can_TP53_9  .   GENE_ID=TP53
chr17   7579853 7579943 Can_TP53_10 .   GENE_ID=TP53
chr2    29445209    29445298    Can_ALK_22  .   GENE_ID=ALK
chr2    29443574    29443664    Can_ALK_23  .   GENE_ID=ALK
chr2    29432659    29432741    Can_ALK_25  .   GENE_ID=ALK
```

### 2、流程跑完，结果目录过大，查看原因是中间的bam文件保留了两次。因此删除*.clean.bam和*tmap.bam，保留过滤比对排序后并用于call变异的*_rawlib.bam。

### 3、Can_28插件修复问题如下，一直显示Started，查看log文件，最后的SGE状态为25。

![](/media/202301/2023-01-09_140820&163vXYMpZdaeoEPfD75x.png)

![](/media/202301/2023-01-09_140813&4SeNADnZ32kgasHExbB9.png)

查找原因：在目录/path/to/Can_28_out.553/BamStat/mpileup下有can28.sh、can28.log和filter.sh、filter.log，在filter.log文件中查看到如下报错：
```
> 报错1
Argument "1stVariationMAF" isn't numeric in numeric gt (>) at /results/plugins/Can_28/myscript/jwli/filter3.pl line 51, <IN> line 1.
#报错2
Can't use an undefined value as an ARRAY reference at /results/plugins/Can_28/myscript/jwli/anno_mpileup.pl line 87, <IN> line 4.
```

报错1在filter3.pl文件中，因为有表头，因此添加内容: `next if /^#Chromosome/; #add by lhlv 20180903`

报错2在anno_mpileup.pl，因此拷贝并修改该文件anno_mpileup2.pl，增加了结果输出文件Tag**-Tag**.3.anno2。

最后，在/results/plugins/Can_28/myscript/jwli/filtersh.py文件中调用相应的脚本，下面为修改后生成的filter.sh文件
```
/results/plugins/Can_28/myscript/jwli/filter3.pl /results/analysis/output/Home/Auto_user_2456177-5-0049-108-P5-HiQ-pooling-Can28-BRCA-CHMO--20180828_240_236/plugin_out/Can_28_out.553/BamStat/mpileup/Tag117.mppana.txt /results/analysis/output/Home/Auto_user_2456177-5-0049-108-P5-HiQ-pooling-Can28-BRCA-CHMO--20180828_240_236/plugin_out/Can_28_out.553/BamStat/mpileup/Tag115.mppana.txt > /results/analysis/output/Home/Auto_user_2456177-5-0049-108-P5-HiQ-pooling-Can28-BRCA-CHMO--20180828_240_236/plugin_out/Can_28_out.553/BamStat/mpileup/04.variant/Tag115-Tag117.3
/results/plugins/Can_28/myscript/jwli/anno_mpileup2.pl /results/analysis/output/Home/Auto_user_2456177-5-0049-108-P5-HiQ-pooling-Can28-BRCA-CHMO--20180828_240_236/plugin_out/Can_28_out.553/BamStat/mpileup/04.variant/Tag115-Tag117.3
```

### 4、bug
![](/media/202301/2023-01-09_140937&MtQsdma24SPXoufGCech.png)

从头一步一步寻找bug，最终发现是在文件uploder/can28_mutation.txt中的里有exonic;splicing这样的变异注释，对应到后一列为两个基因名称了（如TP53;TP53），因为ANNOVAR注释到两种不同的转录本上。而程序中对这样的变异直接读取gene那一列，到最后将分号替换为tab时将会有两列基因出现，从而导致列的错误，从而提取不到相应的变异类型（如nonsysmous），导致hash报错。
因为我们最后只要一个转录本的结果输出，而且发现了这个流程并没有对splicing剪切位点进行注释（ANNOVAR版本过旧，2013年），所以对该bug只要在drug_annotation.py文件中，获取gene时，不直接获取基因那一列，而是通过函数get_all_mutation_type直接拆分AAChange.refgene列的信息中获取gene。

附：2013版ANNOVAR注释结果：（没有GeneDetail.refGene那一列）
![](/media/202301/2023-01-09_140942&2o5YBa1WkIXJtN0vpjZl.png)

最新版（2017年）ANNOVAR注释结果（摘自单病插件）：（GeneDetail.refGene列是标书内含子上信息的）
![](/media/202301/2023-01-09_140947&cgsNF6y3z2bMQmhD7fAI.png)