搜索
查看: 1247|回复: 12

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

[复制链接]

631

主题

1159

帖子

3884

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3884
发表于 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个文件,我

如果你觉得下载文件比较麻烦,我给你一个测试数据如下:
[AppleScript] 纯文本查看 复制代码
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

给出perl代码如下:
[Perl] 纯文本查看 复制代码
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;




上一篇:为什么总是要我修改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
回复 支持 3 反对 0

使用道具 举报

0

主题

1

帖子

129

积分

注册会员

Rank: 2

积分
129
发表于 2017-5-3 21:09:02 | 显示全部楼层
Python
用的是健明老师给的测试数据
用到的方法都是dongye以前讲过的,自己也上网查了一些
[Python] 纯文本查看 复制代码
#!/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
最后将各个染色体的情况分文件输出
手动进行了检查,代码可以使用
如果大家发现我的代码有什么可以提高效率的方法,请告诉我,十分感谢!

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

使用道具 举报

64

主题

137

帖子

657

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
657
发表于 2017-4-28 16:11:06 | 显示全部楼层
R
用的是群主给的测试数据
1.将22条染色体的名称变为一个向量
2.设置循环,用which函数提取某条染色体的所有的行
3.写文件,命名为改染色体名称
[Python] 纯文本查看 复制代码
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")
}

结果:
生成22个文件,例如名为chr7的文件例就是对应的3行chr7的数据。

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

5

帖子

99

积分

注册会员

Rank: 2

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

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

使用道具 举报

25

主题

100

帖子

760

积分

高级会员

Rank: 4

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

见识到了use FileHandle的用法

[Perl] 纯文本查看 复制代码
#!/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;

回复 支持 反对

使用道具 举报

0

主题

16

帖子

145

积分

注册会员

Rank: 2

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


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



回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-5-22 17:50:31 | 显示全部楼层
用群主给的例子数据:
[Python] 纯文本查看 复制代码
#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')


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

使用道具 举报

0

主题

15

帖子

147

积分

注册会员

Rank: 2

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



代码基本上就是原样,如下:
[Perl] 纯文本查看 复制代码
#!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;




本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

8

帖子

117

积分

注册会员

Rank: 2

积分
117
发表于 2017-6-27 19:25:20 | 显示全部楼层
[Python] 纯文本查看 复制代码
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)
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-10-20 03:57 , Processed in 0.135351 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.