搜索
查看: 1744|回复: 0

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

[复制链接]

103

主题

133

帖子

854

积分

版主

Rank: 7Rank: 7Rank: 7

积分
854
发表于 2017-3-10 21:56:54 | 显示全部楼层 |阅读模式
【直播】我的基因组61:scalpel软件找indel
昨天,小编不认真,把微信直播记错了数,感谢认真的小伙伴的发现提醒,你这么认真,一定很快就步入大神的行列的。

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


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

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


Scalpel is available here: [url=]http://scalpel.sourceforge.net/[/url]
文章:
[url=]http://www.nature.com/nmeth/journal/v11/n10/full/nmeth.3069.html[/url]
软件说明书写的也比较详细:[url=]http://scalpel.sourceforge.net/manual.html[/url]

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

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


软件安装是:
  • [AppleScript] 纯文本查看 复制代码
    ## 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软件代码如下:

  • [AppleScript] 纯文本查看 复制代码
    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


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

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

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

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

[url=]https://isogg.org/tree/index.html[/url]
[url=]https://www.yfull.com/tree/[/url]
[url=]http://www.ybrowse.org/gbrowse2/gff/[/url]
[url=]http://www.phylotree.org/tree/index.htm[/url]
[url=]https://dna.jameslick.com/mthap/[/url]   

文:Jimmy

图文编辑:吃瓜群众






上一篇:【直播】我的基因组60:CNV初步探索
下一篇:第七章 直线相关与回归分析(生物统计学基础)
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-23 07:24 , Processed in 0.034226 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.