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

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

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

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

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模式。

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

软件安装是:

## 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 这个文件,我就不多说了,前面我们已经讲过了!

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软件代码如下:

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

https://isogg.org/tree/index.html

https://www.yfull.com/tree/

http://www.ybrowse.org/gbrowse2/gff/

http://www.phylotree.org/tree/index.htm

https://dna.jameslick.com/mthap/

文:Jimmy

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-03-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏java一日一条

在Java 8下更好地利用枚举

在我们的云使用分析API中,返回了格式化过的分析数据(这里指生成分析图)。最近,我们添加了一个特性,允许用户选择时间段(最开始只可以按天选择)。问题是,代码中每...

931
来自专栏java一日一条

有经验的Java开发者和架构师容易犯的10个错误(上)

首先允许我们问一个严肃的问题?为什么Java初学者能够方便的从网上找到相对应的开发建议呢?每当我去网上搜索想要的建议的时候,我总是能发现一 大堆是关于基本入门的...

612
来自专栏Jimoer

Java设计模式学习记录-责任链模式

 已经把五个创建型设计模式和七个结构型设计模式介绍完了,从这篇开始要介绍行为型设计模式了,第一个要介绍的行为型设计模式就是责任链模式(又称职责链模式)。

1092
来自专栏纯洁的微笑

看程序员怎么解决食堂排队问题

在学校的时候,我不爱去食堂成功,一是由于暗黑料理,更重要的一点是人太多了,队伍往往从窗口排到了门口,点菜、计算价格、付款三种业务由打饭阿姨一人完成,思维切换忙碌...

661
来自专栏Data Analysis & Viz

手把手教你完成一个数据科学小项目(2):数据提取、IP查询

本系列将全面涉及本项目从爬虫、数据提取与准备、数据异常发现与清洗、分析与可视化等细节,并将代码统一开源在GitHub:DesertsX/gulius-proje...

841
来自专栏IMWeb前端团队

什么鬼,又不知道怎么命名class了

相信写css的人都会遇到下面的问题: 糟糕,怎么命名这个class,好像不太贴切,要是冲突了怎么办,要不要设计成通用一点... 而改别人css代码的时候则会一直...

2298
来自专栏生信小驿站

TCGA生存分析③

TCGA 癌症基因组图谱(TCGA)是国家癌症研究所(NCI)和国家人类基因组研究所(NHGRI)之间的合作,收集了33种癌症类型的大量临床和基因组数据。 整...

1175
来自专栏狮乐园

【译】Understanding SOLID Principles - Dependency Inversion

当我们在读书,或者在和一些别的开发者聊天的时候,可能会谈及或者听到术语SOILD。在这些讨论中,一些人会提及它的重要性,以及一个理想中的系统,应当包含它所包含的...

823
来自专栏编程

使用JavaScript开发一个自修改代码

话说在25年前,我刚刚开始从事软件开发。在工作中,我遇到一个叫Dave的朋友,他曾在一家大型保险公司工作过几年,他的工作重点是开发支持一个名为“个人人寿保险”的...

2517
来自专栏Java学习网

最佳编码实践:搞砸代码的10种方法

 这是一篇提供有效、实用编程方法的程序箴言,作者Susan Harkins是世界最大的技术期刊出版社的主编,具有多年的实践经验;在这篇文章里她重申“最佳编码实践...

2614

扫码关注云+社区