搜索
楼主: Jimmy

生信编程直播第四题:多个同样的行列式文件合并起来

  [复制链接]

0

主题

4

帖子

97

积分

注册会员

Rank: 2

积分
97
发表于 2017-2-15 16:15:00 | 显示全部楼层
本帖最后由 youzicha 于 2017-2-15 16:18 编辑

编号099
首先我的思路是如果仅仅是合并两个文件要怎么操作,然后再在两个文件处理的基础上思考如何将多个文件进行处理,简而言之就是由少到多的思路。
两个文件进行合并很好处理,一般很容易get到思路,一开始我的思路是用循环嵌套,后来我在运行的时候发现如果用这种方法很耗时,效率不高,在看到之前的小伙伴写的代码之后,我突然意识到哈希可能是一个更好的手段,于是我写成了这样:
[Perl] 纯文本查看 复制代码
#! /usr/bin/perl 
use warnings;
use strict;
my ($line,$header,@info,%name,$head,$id);
$head="gene_ID";
open (NAME,"<name.txt")||die "$!";
while ($id=<NAME>){
        chomp $id;
        $head.="\t$id";
        `gunzip $id.txt.gz`;
        open (FILE,"<$id.txt")||die "$!";
        while ($line=<FILE>){
        chomp $line;
        ($header,@info)=split(/\t/,$line);
        $name{$header}.="@info\t";
        }
}
open (FIN,">final_combination.txt")||die "$!";
print FIN "$head\n";
foreach my $i(sort keys %name){
        print FIN "$i\t$name{$i}\n";
}
close NAME;
close FIN;
最后结果如图:

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-2-15 16:37:20 | 显示全部楼层
本帖最后由 y461650833y 于 2017-2-15 16:55 编辑

今天学习了下用R处理本题,用的是随机生成的数据!
[AppleScript] 纯文本查看 复制代码
awk '{print FILENAME"\t"$0}' * >tmp.txt
setwd("E:/perl/test")
a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)
row.names(fpkm)<-fpkm$V2
fpkm<-fpkm[,-1]
names(fpkm)
 [1] "a.txt" "b.txt" "c.txt" "d.txt" "e.txt" "f.txt" "g.txt" "h.txt" "i.txt" "j.txt" "k.txt"
[12] "l.txt" "m.txt" "n.txt" "o.txt" "p.txt" "q.txt" "r.txt" "s.txt" "t.txt" "u.txt" "v.txt"
[23] "w.txt" "x.txt" "y.txt" "z.txt"
names(fpkm)<-unlist(lapply(names(fpkm),function(x){
+   fpkm<-strsplit(x,'\\.')[[1]][1]
+ }))
write.csv(fpkm,"test.csv")



同理,处理GSE42813 也得到对应的结果

027-R+Perl-游戏




本帖子中包含更多资源

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

x
人生若只如初见!
回复 支持 反对

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-15 21:46:42 | 显示全部楼层
002+python
1. 虽然程序很短,但是每一条都是精华。
2. 开始了解输出type() 的重要性。
3. 还有一些不输入的默认参数需要牢记,例如append()默认加在最后一个,readline(),默认一次读一行。
4. 字典的键值只能有一个,但后面的参数可以有很多。



[Python] 纯文本查看 复制代码
#! /usr/bin/python

import glob  # 查找文件系统中指定模式的路径

MyDiction = {}

def Readfile(File):    #这个函数面向的对象是FILE
    file_obj = open(File)
    try:              #try...finally 记住这个格式!
        while True:
            line = file_obj.readline().strip("\n")   #返回一个str,()默认一次读一行
            if not line:
                break
            array = line.split("\t")      #返回一个两个变量的list
            if array[0] in MyDiction:
                MyDiction[array[0]].append(array[1])  # append()默认加在后面
            else:
                MyDiction[array[0]] = [array[1]]  #如果没有,第一次录入键+参数
    finally:
        file_obj.close()
    return MyDiction

def main():
    list_dirs = glob.glob("./*.txt")   #glob 包里的glob 函数 引入当前目录下.txt 结尾的文件
    for i in list_dirs:
        Readfile(i)
    for gene_name in MyDiction:
        print ("%s\t%s" %(gene_name, "\t".join(MyDiction[gene_name])))    

if __name__ == '__main__':   # 文件运行的入口,从主函数开始
    main()


基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 3 反对 0

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-15 21:53:56 | 显示全部楼层
熊海清 发表于 2017-2-6 22:10
074-python
用的是论坛帖子中做WGCNA分析的56个sample,主要是理解利用字典存储结果,且将数组作为字典的值 ...

我水平和你差不多,但我感觉你理解可能有误差。
line = file_obj.readline().strip("\n")        这个返回的是去掉/n  的字符串     gene_1 123/n →    gene_1 123
array = line.split("\t")                          这返回一个list                         gene_1 123    →  ['gene_1','123']
Mydic[array[0]] = [array[1]]                这是一个key 键不是赋值      ['gene_1','123']   → {gene_1:123}
我是新手,欢迎交流
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

0

主题

2

帖子

40

积分

新手上路

Rank: 1

积分
40
发表于 2017-2-16 22:31:27 | 显示全部楼层
从 这次的作业中学到了很多东西,尤其是读入多个文件的方式;还有就是join函数的用法,这些是之前自己没有关注到的点,学习就是个积累的过程,加油加油!
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my %hash;
my $header="Gene_name";
foreach my $i(0..@ARGV-1) {   ##从这行 学到了一种读入多个文件的方式
        $ARGV[$i]=~/(.*?)\.txt/;
        $header.="\t$1";
        open IN,"$ARGV[$i]";
        while  (<IN>)   {
                chomp;
                my @t=split;
                $hash{$t[0]}[$i]=$t[1];
        }
}
print $header,"\n";
foreach my $k (sort keys %hash)  {
        foreach my $i (0..@ARGV-1)  {
                $hash{$k}[$i] =0  if !defined($hash{$k}[$i]);
        }
        print  $k,"\t",join("\t",@{$hash{$k}}),"\n";
}        
附下部分结果:
Gene_name       GSM1172844_184A1        GSM1172845_184B5        GSM1172846_21MT1        GSM1172847_21MT2        GSM1172848_21NT GSM1172849_21PT GSM1172850_600MPE       GSM1
ENSG00000000003 95.2125479443   95.6986763046   19.9946736268   65.6863762916   44.0577455869   34.3175652587   178.1588315834  13.4601437726   37.8796664069   69.092581378
ENSG00000000005 0.0000000000    0.0000000000    0.0000000000    0.1492020573    0.0000000000    0.0000000000    0.0000000000    0.0000000000    0.0000000000    0.0000000000
ENSG00000000419 453.2083073388  243.6480387207  142.0581753071  200.4131492951  193.1543893082  151.5729124615  220.7534899930  147.4689279686  684.5909386275  314.03122029
ENSG00000000457 18.1043900660   26.5666075008   16.1277638108   12.0873135093   18.4480093004   15.7835259267   89.6726759246   34.4605135027   140.4273980447  70.603729347
ENSG00000000460 48.1662237375   24.5842890228   24.2845922387   36.5169168453   32.5867631826   28.5225518204   140.5434281494  16.7714033080   92.5539077330   46.691517271
回复 支持 反对

使用道具 举报

0

主题

6

帖子

89

积分

注册会员

Rank: 2

积分
89
发表于 2017-2-19 20:56:18 | 显示全部楼层
本帖最后由 banapi 于 2017-2-21 01:54 编辑

158
shell
shell中使用perl,perl中","是运算符,"."是连接符
-M5.010能增加一些语法,不加这个的话不能用say,say会自动加换行
paste只能简单的列合并
join可以关联合并,但只能文件两两进行

构建文件集:
[Shell] 纯文本查看 复制代码
perl -M5.010 -e 'foreach(a..z){say $_."\.txt"}'|while read line; do perl -M5.010 -e 'foreach(1..99){say "gene_$_\t".int(rand(1000))}'>$line; done

合并:
[Shell] 纯文本查看 复制代码
perl -M5.010 -lane 'BEGIN{$a} push @{$a->{$F[0]}},$F[1];END{foreach(sort keys %$a){say $_,"\t",(join "\t",@{$a->{$_}})}}' *txt|less -S

部分结果
gene_1  622     926     18      3       89      627     76      116     325     571     237     744     455     285     693     551
gene_10 900     425     467     321     918     31      208     820     682     500     849     549     691     946     801     95
gene_11 282     758     387     816     907     473     305     236     845     201     724     949     620     950     629     352
gene_12 338     486     667     834     313     580     252     229     483     757     622     379     163     111     483     781
gene_13 699     2       202     87      385     28      194     887     267     471     537     630     824     581     640     501
gene_14 222     213     195     671     64      366     479     534     677     272     555     48      651     344     793     732
gene_15 439     296     834     169     202     184     565     759     508     323     23      80      50      565     481     120
gene_16 777     164     720     156     229     82      241     677     667     284     181     94      879     318     161     224








回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-2-19 23:18:43 | 显示全部楼层
本帖最后由 xxnku 于 2017-5-7 21:17 编辑

biotrainee_four.py
[Python] 纯文本查看 复制代码
import glob
from collections import OrderedDict

tmp_dict = OrderedDict()

def readfile(File):
    try:
        with open(File) as f:
            for line in f:
                line = line.rstrip()
                if not line:
                    continue
                gene_name = line.split('\t')[0]
                gene_num = line.split('\t')[1]
                if gene_name not in tmp_dict.keys():
                    tmp_dict[gene_name] = [gene_num]
                else:
                    tmp_dict[gene_name].append(gene_num)
    finally:
        pass
    return tmp_dict

def main():
    list_dir = glob.glob('./*.txt')

    for i in list_dir:
        readfile(i)

    for k,v in tmp_dict.items():
        print ('%s\t%s'%(k,'\t'.join(v)))

if __name__ == '__main__':
    main()

python3 biotrainee_four.py | less -S
回复 支持 反对

使用道具 举报

0

主题

8

帖子

179

积分

注册会员

Rank: 2

积分
179
发表于 2017-2-20 23:57:13 | 显示全部楼层
028-Python-Ryu

这段时间一直忙其他的东西和偷懒,觉得很久不练习果然落后很多,第4题其实很简单于是赶紧捡起来完成。

一开始以为这一题的每个文件中的基因是不一样的,有的有,而有的没有,就觉得会复杂很多,也是我看题不仔细;后来看了老师的视频后,发现原来都是对齐的。

解题思路:
1.使用python的gzip模块读取tar解压后的56个文件,减少空间以及解压的时间;
2.使用python的os模块中的os.walk遍历所有文件,在遍历中写入字典;
3.使用字典存储每个基因的每个文件中的数值,使用字典类似哈希,使用键值对方式存储,复杂度低,索引速度快;
4.使用Python的collection模块中的defaultdict,该字典支持字典嵌套,可以很方便地写出嵌套的字典,容易理解;
5.使用Python中的csv模块读写文件,使用dictReader读文件,同文件中每行的key相同,同时将每个文件中的列标题存入一个列表中,因为Python中的字典是无序的,使用列表来索引可以保持有序,同时用作输出文件的标题
6.使用Win系统时,要注意gzip.open读取文件要用‘rt’,csv.writer要在打开文件时设置'newline=''‘,原因之前某同学说过,CRLF与LF的问题。

代码如下:
[Python] 纯文本查看 复制代码
#!/usr/bin/python
# -*- coding: utf-8 -*-
'''
/*
 * @Author: Ryu.Zheng 
 * @Date: 2017-02-19 17:55:49 
 * @Last Modified by:   Ryu.Zheng 
 * @Last Modified time: 2017-02-19 17:55:49 
 */
'''

import os
from collections import defaultdict
import gzip
import csv

gene_table=defaultdict(dict)
##NOTE: {gene_table: 
#           {Gene_ID:{filename:value},{filename:value},...},
#           {Gene_ID:{filename:value},{filename:value},...},
#           ...
#           }

file_list=[]

for root,dirs,files in os.walk("GSE48213_RAW"):
    for name in files:
        # print(name)
        with gzip.open(os.path.join("GSE48213_RAW",name),'rt') as gene_file:
            header=next(csv.reader(gene_file,delimiter="\t"))
            print(header)
            file_list.append(header[1])

            genes = csv.DictReader(gene_file,delimiter="\t",fieldnames=header)
            for line in genes:
                # print(line)
                if line[header[0]] in gene_table:
                    gene_table[line[header[0]]][header[1]]=line[header[1]]
                else:
                    gene_table[line[header[0]]]={header[1]:line[header[1]]}


with open("merge_table.txt",'w',newline="")as f2:
    csvwrite=csv.writer(f2,delimiter="\t")
    csvwrite.writerow(["Genes"]+file_list)

    for gene in sorted(gene_table):
        ## header
        # for name in gene_table[gene]:
        content=[gene_table[gene][name] for name in file_list]
        csvwrite.writerow([gene]+content)


欢迎各位同学交流,谢谢老师们辛苦的劳动,我会尽快补齐作业,多写博客交流学习心得
回复 支持 反对

使用道具 举报

0

主题

5

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-2-24 18:08:26 | 显示全部楼层
本帖最后由 chzhniu 于 2017-2-24 18:18 编辑

182-perl-牛

首先是下载文件并对文件进行解压,学习Linux文件压缩与解压不同方法ta、gz等。
以下是根据授课内容然后按照自己的理解写的代码:
[Perl] 纯文本查看 复制代码
#! user/bin/perl
use strict;

my %total;
my $header = "gene";
foreach my $i(0..@ARGV-1){ #@ARGV是由命令参数组成的数组,即要合并的文件
$ARGV[$i]=~ /(.*?).txt/; #对文件名进行匹配
$header .="\t$1"; #将匹配得到的文件名作为合并后文件的文件的表头

open I, "$ARGV[$i]";
while(<I>){
        chomp;
        my @tmp = split; 
        my $k = $tmp[0];
        $total{$k}[$i] = $tmp[1]; #对每一个文件进行拆分,第一列为%total哈希的键,第二列为对应的值,并转化为数组@{$total{key}}
        }
}
print "$header\n"; #打印表头
close I;

foreach my $key (sort keys %total){
        foreach my $i(0..@ARGV-1){ #该处循环的目的是看看不同文件中的键值(对应的基因)是否相同,如果有缺少的,将缺省值赋值为0,如果没有这个循环,结果中有缺省的值会是空的
        $total{$key}[$i]= 0 if !defined $total{$key}[$i];
        }
print $key, "\t", join ("\t", @{$total{$key}}), "\n"; #打印结果$key表示所有基因,join将不同数组结合起来
}


通过本题主要学习合并多个同样行列文件的方法,同时学习多维数组哈希的用法(还需进一步理解学习)。
回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-2-26 18:04:41 | 显示全部楼层
zckoo007 发表于 2017-2-15 21:46
002+python
1. 虽然程序很短,但是每一条都是精华。
2. 开始了解输出type() 的重要性。

感谢你呀,你的注释让我看懂了好多!
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-13 03:56 , Processed in 0.038224 second(s), 20 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.