scalpel软件找indel

Posted on by

Scalpel is available here: http://scalpel.sourceforge.net/
文章是： http://www.nature.com/nmeth/journal/v11/n10/full/nmeth.3069.html
很赞的工具！
软件说明书写的也比较详细：http://scalpel.sourceforge.net/manual.html
他提供了3种情况的找INDELs变异，我目前需要用的就是对我的全基因组测序数据来找，所以用single模式：
为了节省对计算资源的消耗，作者建议我单独对每条染色体分别处理。

软件安装是：

1
2
3
4
5
6
7
8
9
## Download and install Scalpel
cd ~/biosoft
mkdir Scalpel &&  cd Scalpel
wget [url]https://downloads.sourceforge.net/project/scalpel/scalpel-0.5.3.tar.gz[/url] 
tar zxvf scalpel-0.5.3.tar.gz
cd scalpel-0.5.3
make
~/biosoft/Scalpel/scalpel-0.5.3/scalpel-discovery  --help
~/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软件代码如下：

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
cat ~/reference/genome/hg19/hg19.chr.bed |while read id
do
arr=($id)
# arr=($a) will split the $a to $arr , ${arr[0]} ${arr[1]} ~~~, but ${arr[@]}  is the whole array .
# OLD_IFS="$IFS"
# IFS=","
# arr=($a)
# IFS="$OLD_IFS"
#arr=($a)用于将字符串$a分割到数组$arr ${arr[0]} ${arr[1]} ... 分别存储分割后的数组第1 2 ... 项 ，${arr[@]}存储整个数组。
#变量$IFS存储着分隔符，这里我们将其设为逗号 "," OLD_IFS用于备份默认的分隔符，使用完后将之恢复默认。
echo ${arr[0]}:${arr[1]}-${arr[2]}
 
date
start=`date +%s`
~/biosoft/Scalpel/scalpel-0.5.3/scalpel-discovery --single \
--bam  ~/data/project/myGenome/fastq/bamFiles/jmzeng.filter.rmdup.bam \
--ref ~/reference/genome/hg19/hg19.fa \
--bed ${arr[0]}:${arr[1]}-${arr[2]}  \
--window 600 --numprocs 5  --dir ${arr[0]}
end=`date +%s`
runtime=$((end-start))
echo "Runtime for ${arr[0]}:${arr[1]}-${arr[2]} was $runtime"
done

Comments are closed.