搜索
查看: 8064|回复: 16

生信编程直播第11题:把文件内容按照染色体分开写出

[复制链接]

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-4-11 16:48:32 | 显示全部楼层 |阅读模式
这个需求很常见,因为一般生物信息学数据比较大,比如sam,vcf,或者gtf,bed都是把所有染色体综合在一起的文件。
如果想根据染色体把大文件拆分成小的文件呢?
比如:ftp://ftp.ncbi.nlm.nih.gov/pub/C ... an/CCDS.current.txt
这个文件里面的基因就有染色体信息,那根据染色体把这个文件拆分成1~22和其它染色体,这样的23个文件,我

如果你觉得下载文件比较麻烦,我给你一个测试数据如下:
[mw_shl_code=applescript,true]chr2        43995310        43995986
chr17        49788603        49789067
chr17        59565573        59566163
chr19        8390308        8390745
chr12        49188033        49189033
chr7        974903        975570
chr7        98878532        98879500
chr7        44044672        44045322
chr1        153634052        153634772
chr11        60905850        60906575
[/mw_shl_code]
给出perl代码如下:
[mw_shl_code=perl,true]use FileHandle;
foreach( 1..22 ){
        $normal_chr{"chr".$_}=1 ;
        open $fh{"chr".$_},">chr$_.txt" or die;
}
open other,">chr_other.txt" or die;
open FH,'a.bed';
while(<FH>){
        chomp;
        @F=split;
        if(exists $normal_chr{$F[0]}){
                 $fh{$F[0]}->print( "$_\n" );
        }else{
                print other $_;
        }
}
foreach( 1..22 ){
        close $fh{$_};
}
close other;[/mw_shl_code]



上一篇:为什么总是要我修改mysql密码?
下一篇:从metadata中提取临床信息
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

0

主题

2

帖子

35

积分

新手上路

Rank: 1

积分
35
发表于 2017-4-12 09:32:49 | 显示全部楼层
awk '{print >> ("CCDS.current.txt"$1)}’ CCDS.current.txt
回复 支持 4 反对 0

使用道具 举报

0

主题

1

帖子

139

积分

注册会员

Rank: 2

积分
139
发表于 2017-5-3 21:09:02 | 显示全部楼层
Python
用的是健明老师给的测试数据
用到的方法都是dongye以前讲过的,自己也上网查了一些[mw_shl_code=python,true]#!/usr/bin/env python

import sys

args = sys.argv

targ = args[1]

with open(targ, 'r') as F:
        Chr = {}
        for line in F:
                line = line.split('\t', 1)
                if line[0] not in Chr:
                        Chr[line[0]] = ''
                        Chr[line[0]] += line[1]
                else:
                        Chr[line[0]] += line[1]

list = Chr.keys()
n = len(list)

i = 1

while i < n+1:
        with open('out%s.txt'%(i), 'w') as out:
                out.write('%s\t'%(list[i-1]))
                out.write('%s'%(Chr[list[i-1]]))

        i += 1[/mw_shl_code]最后将各个染色体的情况分文件输出
手动进行了检查,代码可以使用
如果大家发现我的代码有什么可以提高效率的方法,请告诉我,十分感谢!

在作业刚刚出来的时候自己就已经想到了怎么写,也很快就写完了;
今天突然觉得自己应该分享出来,我这个学期才进入生信菜鸟团,之前也只是学了Python的基础语句,这么多编程实战,我从第一个开始,一直都是看着视频,听老师讲一遍,然后自己一点一点敲,在自己的电脑上实现。这次作业的时候,我发现自己竟然可以独立做出来了,虽然仍然是个菜鸟,但自己还是有点小激动的。
感谢Jimmy和各位不吝啬拿出自己时间来分享经验的老师,真的学到很多。
回复 支持 1 反对 0

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-4-28 16:11:06 | 显示全部楼层
R
用的是群主给的测试数据
1.将22条染色体的名称变为一个向量
2.设置循环,用which函数提取某条染色体的所有的行
3.写文件,命名为改染色体名称
[mw_shl_code=python,true]rm(list=ls)
install.packages("stringr")
library("stringr")
setwd("D:\\myR")
a<-read.table(file="test.txt",header = FALSE)
name<-c("chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8","chr9","chr10","chr11",
        "chr12","chr13","chr14","chr15","chr16","chr17","chr18","chr19","chr20","chr21","chr22")
for(n in 1:22){
b<-a[which(a[,1]%in%c(name[n])),]
write.table(b,file=name[n],row.names=FALSE,sep="\t")
}[/mw_shl_code]
结果:
生成22个文件,例如名为chr7的文件例就是对应的3行chr7的数据。

无标题.png
回复 支持 反对

使用道具 举报

0

主题

5

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-4-29 22:33:57 | 显示全部楼层
本帖最后由 王土土哈 于 2017-4-30 14:22 编辑

[mw_shl_code=shell,true]for i in `cut -f1 CCDS.current.txt|sort -k1,1 -u|grep -v "#"`; do awk -v id=$i '{if($1==id)print $0}' CCDS.current.txt >Chr${i}.txt; done
awk '($1!~"#"){print > "Chr"$1"."FILENAME}' CCDS.current.txt[/mw_shl_code]
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-5-4 22:40:52 | 显示全部楼层
034-perl-shell

见识到了use FileHandle的用法

[mw_shl_code=perl,true]#!/usr/bin/perl -w
use strict;
use FileHandle;

my %hash;
my %tmp;
map{
        $hash{$_} =1;
        open $tmp{$_}, ">chr$_.txt" or die "Cannot write it";
}(1..20);

open my $fh_other, ">chr_other.txt" or die "Cannot write chr_other.txt";
open my $fh, "CCDS.current.txt" or die "Cannot open CCDS.current.txt";
readline $fh;
while(<$fh>){
        chomp;
        my @array = split /\t/, $_;
        if($hash{$array[0]}){
                print "aaa\n";
                $tmp{$array[0]} -> print ("$_\n");
        }else{
                print $fh_other ("$_\n");
        }
}

map{
        close $tmp{$_};
}(1..20);
close $fh_other;
[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

16

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-15 19:59:22 | 显示全部楼层
相对来说比较简单,直接调用第零题的readfasta函数。os包不是很熟悉,只能现学现用了。
觉得输出一堆文件看着麻烦,就没有实际运行。


[mw_shl_code=applescript,true]def readfasta(filename):
    """
    :param filename: ;要读取的FASTA文件名
    :return: 返回基因名:序列的字典文件
    """

    fa = open(filename, 'r')
    res = {}
    ID = ''
    for line in fa:
        if line.startswith('<'):
            ID = line.strip('\n')
            res[ID] = ''
        else:
            res[ID] += line.strip('\n')

    return res

def seperatefile(filename, outfile):
    data=readfasta(filename)
    for key, value in data.items():
        import os
        (file, ext)=os.path.splitext(outfile)
        fo=open('%s_%s%s' % (file, key, ext), 'w')
        fo.write('%s\n'%key)
        fo.write('%s\n'%value)


seperatefile('testfa.fa', 'output.fa')
[/mw_shl_code]


回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-5-22 17:50:31 | 显示全部楼层
用群主给的例子数据:
[mw_shl_code=python,true]#encoding = utf-8

from collections import OrderedDict

def readfasta(filename):

    tmp_dict = OrderedDict()

    with open(filename) as f:

        for line in f:
            line = line.rstrip().split(' ',1)

            chr_id = line[0]

            if chr_id not in tmp_dict:
                tmp_dict[chr_id] = line[1]

            else:
                tmp_dict[chr_id] += line[1]

    return tmp_dict

def seperatefile(filename,outfile):

    data = readfasta(filename)

    for chr_id,features in data.items():
        import os
        (name,ext) = os.path.splitext(outfile)

        with open('%s_%s%s' %(name,chr_id,ext),'w') as f_out:
            f_out.write('%s\n' %chr_id)
            f_out.write('%s\n' %features)

seperatefile('test.txt','output.txt')[/mw_shl_code]

谢谢dongye和azazelcc,学到了os.path
回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-6-8 15:07:16 | 显示全部楼层
做这道题对于我来说,感觉最大的难点在于怎么进行批量文件生成,直到仔细看了楼主的代码,我真是被惊吓掉了下巴,原来除了解引用还有一个类的用法,学习了,这个方法我这连天在持续联系,感觉很好用啊;还有就是,use FileHandle被我注释掉为什么还是可以正常工作,这一点没有太明白



代码基本上就是原样,如下:
[mw_shl_code=perl,true]#!usr/bin/perl
use warnings;
use strict;
#use FileHandle;

my (%chr,%fh);

foreach( 1..22 ){
        $chr{$_}=1 ;
        open $fh{$_},">chr$_.txt" or die;
}

open Other,">chr_other.txt" or die;
open FH,$ARGV[0] or die $!;

while(<FH>){
        chomp;
        my @info=split;
        if(exists $chr{$info[0]}){
                $fh{$info[0]}->print( "chr$_\n" );#这个用法我真的涨姿势了;
        }else{
                print Other "$_\n";
        }
}

foreach( 1..22 ){
        close $fh{$_};
}

close FH;
close Other;[/mw_shl_code]



运行信息

运行信息

生成文件

生成文件

Chr11的结果展示

Chr11的结果展示
回复 支持 反对

使用道具 举报

0

主题

8

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-6-27 19:25:20 | 显示全部楼层
[mw_shl_code=python,true]import sys
args=sys.argv



def chrseq(inputfile,outputfile):
        with open(inputfile) as F_in:
                with open(outputfile,'w') as F_out:
                        for line in F_in:
                                if line.split('\t')[0]==args[3]:
                                        print(line)
                                        F_out.write(line)
                    
        return outputfile
def main(args):
        inputfile=args[1]
        outputfile=args[2]
        chrseq(inputfile, outputfile)
if __name__=="__main__":
        main(args)[/mw_shl_code]
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-3-28 17:21 , Processed in 0.032558 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.