搜索
查看: 2534|回复: 2

生信编程直播第11题:区分染色体分别运行scalpel软件!(shell)

[复制链接]

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-6 09:17:01 | 显示全部楼层 |阅读模式
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模式:
为了节省对计算资源的消耗,作者建议我单独对每条染色体分别处理。



软件安装是:
[mw_shl_code=applescript,true]## Download and install Scalpel
cd ~/biosoft
mkdir Scalpel &&  cd Scalpel
wget https://downloads.sourceforge.ne ... calpel-0.5.3.tar.gz  
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
[/mw_shl_code]


它需要自己指定--bed参数来选择染色体运行,而且不是给一个chr1就可以了,需要指定染色体及其起始终止坐标:single region in format chr:start-end (example: 1:31656613-31656883)
所以就考验shell编程技巧啦!
制作 ~/reference/genome/hg19/hg19.chr.bed  这个文件,我就不多说了,前面我们已经讲过了!
[mw_shl_code=applescript,true]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
[/mw_shl_code]

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

[mw_shl_code=shell,true]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 [/mw_shl_code]



上一篇:【直播】我的基因组 41:按照不同的lane来call variation
下一篇:3000多份水稻全基因组测序数据
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

3

帖子

58

积分

注册会员

Rank: 2

积分
58
发表于 2017-4-21 19:26:29 | 显示全部楼层
做bed文件过程可否介绍一下
回复 支持 反对

使用道具 举报

7

主题

32

帖子

216

积分

中级会员

Rank: 3Rank: 3

积分
216
发表于 2017-4-28 10:45:12 | 显示全部楼层
本帖最后由 babytong 于 2017-4-28 11:01 编辑

1.python pandas 可以轻松实现文件的提取与输出。
[mw_shl_code=python,true]import os
import pandas as pd
import numpy as np

os.chdir('D:\\')

data=pd.read_table('CCDS.current.txt')
for i in range(1,23):
    df=data.loc[data["#chromosome"]==str(i)]
    df.to_csv('chr'+str(i)+'.csv')
   
df=data.loc[data["#chromosome"]=='X']
df.to_csv('chrX.csv')
df=data.loc[data["#chromosome"]=='Y']
df.to_csv('chrY.csv')[/mw_shl_code]

2.perl 实现此功能也不难
[mw_shl_code=perl,true]use strict;
use warnings;
open(IN,"D:\\CCDS.current.txt") or die "can't open1 data $!\n";
<IN>;
my $col_index=0;  #m代表按照第几列的值进行拆分
my $line=<IN>;
chomp($line);
my $sum=0;
my @a=split(/\s+/,$line);
my $temp=$a[$col_index];
open(OUT,">>D:\\chr_$temp.txt") or die "can't open2 data $!\n";
print OUT "$line\n";
$sum++;
while(<IN>){
        chomp($_);
        my @b=split(/\t/,$_);
        if($b[$col_index] eq $temp){
                print OUT "$_\n";
                $sum++;
        }
        else{
                close(OUT);
                print "$sum\n";
                open(OUT,">>D:\\chr_$b[$col_index].txt")  or die "can't open3 data $!\n";
        print OUT "$_\n";
        $sum++;
                $temp=$b[$col_index];               
                }
                }
        close(IN);
        close(OUT);
        print "$sum\n";[/mw_shl_code]


025 python perl

result_perl.png
resultpython.png
欢迎关注https://www.jianshu.com/u/4f5e357a6212
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-1-20 21:37 , Processed in 0.026003 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.