搜索
查看: 2177|回复: 2

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

[复制链接]

634

主题

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模式:
为了节省对计算资源的消耗,作者建议我单独对每条染色体分别处理。



软件安装是:
[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]https://downloads.sourceforge.ne ... calpel-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  这个文件,我就不多说了,前面我们已经讲过了!
[AppleScript] 纯文本查看 复制代码
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软件代码如下:

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



上一篇:【直播】我的基因组 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

主题

34

帖子

216

积分

中级会员

Rank: 3Rank: 3

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

1.python pandas 可以轻松实现文件的提取与输出。
[Python] 纯文本查看 复制代码
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')


2.perl 实现此功能也不难
[Perl] 纯文本查看 复制代码
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";



025 python perl

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

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

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 09:12 , Processed in 0.033983 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.