【直播】我的基因组61：scalpel软件找indel

昨天，小编不认真，把微信直播记错了数，感谢认真的小伙伴的发现提醒，你这么认真，一定很快就步入大神的行列的。

那么现在正式的开始第61讲：

其实这次的call variation的软件，不仅仅是找到SNV，也顺便找到了indel，只是可能不太准确。一般业界的公认标准是 GATK的best practice，不过那个我已经做了，现在来一点新的，我正好看到了这个scalpel软件。

当然，为什么使用它，完全是随心所欲，也可以选择Pindel等其它软件。我在这里只是为了秀一个软件的用法，生信工程师该如何持续学习。

Scalpel is available here: 

文章：

软件说明书写的也比较详细：

他提供了3种情况的找INDELs变异，我目前需要的就是对我的全基因组测序数据来找，所以用single模式。

为了节省对计算资源的消耗，作者建议我单独对每条染色体分别处理。

软件安装是：

1. ## Download and install Scalpel

2. cd ~/biosoft

3. mkdir Scalpel && cd Scalpel

4. wget [url]https://downloads.sourceforge.net/project/scalpel/scalpel-0.5.3.tar.gz[/url]

5. tar zxvf scalpel-0.5.3.tar.gz

6. cd scalpel-0.5.3

7. make

8. ~/biosoft/Scalpel/scalpel-0.5.3/scalpel-discovery --help

9. ~/biosoft/Scalpel/scalpel-0.5.3/scalpel-export --help

它需要自己指定--bed参数来选择染色体运行，而且不是给一个chr1就可以了，需要指定染色体及其起始终止坐标：single region in format chr:start-end (example: 1:31656613-31656883)，所以就比较考验shell编程技巧啦！

制作 ~/reference/genome/hg19/hg19.chr.bed 这个文件，我就不多说了，前面我们已经讲过了！

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22

chr10 1 135534747 chr11 1 135006516 chr12 1 133851895 chr13 1 115169878 chr14 1 107349540 chr15 1 102531392 chr16 1 90354753 chr17 1 81195210 chr18 1 78077248 chr19 1 59128983 chr1 1 249250621 chr20 1 63025520 chr21 1 48129895 chr22 1 51304566 chr2 1 243199373 chr3 1 198022430 chr4 1 191154276 chr5 1 180915260 chr6 1 171115067 chr7 1 159138663 chr8 1 146364022 chr9 1 141213431

区分染色体分别运行scalpel软件代码如下：

1. cat ~/reference/genome/hg19/hg19.chr.bed |while read id

2. do

3. arr=($id)

4. # arr=($a) will split the $a to $arr , ${arr[0]} ${arr[1]} ~~~, but ${arr[@]} is the whole array .

5. # OLD_IFS="$IFS"

6. # IFS=","

7. # arr=($a)

8. # IFS="$OLD_IFS"

9. #arr=($a)用于将字符串$a分割到数组$arr ${arr[0]} ${arr[1]} ... 分别存储分割后的数组第1 2 ... 项 ，${arr[@]}存储整个数组。

10. #变量$IFS存储着分隔符，这里我们将其设为逗号 "," OLD_IFS用于备份默认的分隔符，使用完后将之恢复默认。

11. echo ${arr[0]}:${arr[1]}-${arr[2]}

12. date

13. start=`date +%s`

14. ~/biosoft/Scalpel/scalpel-0.5.3/scalpel-discovery --single \

15. --bam ~/data/project/myGenome/fastq/bamFiles/jmzeng.filter.rmdup.bam \

16. --ref ~/reference/genome/hg19/hg19.fa \

17. --bed ${arr[0]}:${arr[1]}-${arr[2]} \

18. --window 600 --numprocs 5 --dir ${arr[0]}

19. end=`date +%s`

20. runtime=$((end-start))

21. echo "Runtime for ${arr[0]}:${arr[1]}-${arr[2]} was $runtime"

22. done

最后得到的是每一条染色体一个vcf文件记录着INDEL情况，暂时我还没进行下一步处理。

这里我其实主要是想讲如何用shell进行并行，查看原文可以看到我们的题目及视频讲解，关于这个软件的并行使用！

顺便预告一下，我在wegene测得的芯片数据已经完成了全流程，下载是wegene专题。

还有，我们生信菜鸟团热心群友指出了我前面用常染色体做祖源分析的不足之处，希望我可以继续用Y染色体和线粒体DNA来做下去，给了我几个网址，我估计要学习两个月左右才能完全搞明白，毕竟是孤家寡人兼职学习，有点累，有兴趣的可以学习下面的内容，跟我交流，我的email是[e-mail] 

文：Jimmy

图文编辑：吃瓜群众