搜索
楼主: Jimmy

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

  [复制链接]

0

主题

15

帖子

88

积分

注册会员

Rank: 2

积分
88
发表于 2017-2-11 00:45:02 | 显示全部楼层
#246

完成作业,写得不够lazy啊。。

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

use strict;
use warnings;

opendir DIR, './' or die;
open OUT, '>'.'combined_table';

my (%samples,%genes,%express);

foreach (readdir DIR){
  if(/(.+?)_(.+)\.txt/){
    my $sample = $2;
    $samples{$sample}++;
    open IN, $_ or die;
    while(<IN>){
      chomp;
      next if /Gene_ID/;
      if(/(.+?)\t(.+)/){
        $express{$sample}{$1} = $2;	
        $genes{$1}++;
      }
    }
    close IN;	
  }	
}

my $head = "Gene";
foreach (sort keys%samples){
  $head .= "\t$_";	
}

print OUT "$head\n";

foreach my $gene (sort keys%genes){
  print OUT "$gene";
  my $express;
  foreach my $sample (sort keys%samples){
    $express .= "\t".$express{$sample}{$gene};	
  }
  print OUT "$express\n";
}

回复 支持 反对

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-2-11 08:56:16 | 显示全部楼层
lefroyqiu 发表于 2017-2-11 00:45
#246

完成作业,写得不够lazy啊。。

是的,不够懒,没有得到perl的精髓,但是你的思路是对的,对初学者而言,值得表扬了
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-2-11 10:39:12 | 显示全部楼层
023 R
这次视频看的时间稍晚,其中对字符串的处理和Lapply等函数的应用对我帮助很大
[AppleScript] 纯文本查看 复制代码
#下载安装git,解压56个文件后保存在tmp文件夹下
#启用gitbash 命令awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt 得到汇总文件
setwd("C:\\Users\\qin\\Desktop\\tmp\\tmp")
a<-read.table('tmp.txt')#读入数据
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)#cast是长变宽,v2是行V1是列,默认V3作为value column
rownames(fpkm)<-fpkm[,1]#第一列作为行名
fpkm<-fpkm[,-1]#去掉第一列
names(fpkm)<-unlist(lapply(names(fpkm),function(x){
  tmp<-strsplit(x,'_')[[1]][2]
  strsplit(tmp,'\\.')[[1]][1]
}
))#对字符串的处理,对列名进行精简
write.csv(fpkm," fpkm.csv")
回复 支持 反对

使用道具 举报

0

主题

6

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 2017-2-12 14:44:44 | 显示全部楼层
本帖最后由 雨林木风 于 2017-2-16 11:44 编辑

031-python
思路:
1、遍历每个样本基因表达文件
2、获取每个文件中的基因名,表达量,构建以基因名为keys,不同样本表达量值为value的字典。
3、将字典写入新文件,构建表达矩阵。      
代码:
[Python] 纯文本查看 复制代码
#! usr/bin/env python
# -*- coding:utf-8 -*-

import time
import os
import gzip
import sys
from collections import OrderedDict
start = time.clock()

workdir = sys.argv[1]
os.chdir(workdir) #改变当前工作目录到指定的路径
mergedict = OrderedDict()

#定义函数:读取文件夹内各文件并构建字典用以存储基因及在不同样本的表达量
def readfile(filenames):
    for filename in filenames:
        with gzip.open(filename,'r') as f:
            for line in f:
                line = line.strip()
                line = line.split('\t')
                genename = line[0]
                FPKM = line[1]
                if genename in mergedict:
                    mergedict[genename] += ('\t'+ str(FPKM))
                else:
                    mergedict[genename] = str(FPKM)
    return mergedict

#定义主函数:构建表达矩阵
def main():
    filenames = os.listdir(workdir) #返回指定的文件夹包含的文件或文件夹的名字的列表
    readfile(filenames)
    with open('mergedict.txt','w') as output:
        for genename, FPKM in mergedict.items():
             output.write("%s\t%s\n" %(genename, FPKM))
main()
end = time.clock()
print("usded %s s" % str(end-start))


小结:1、学会了遍历文件夹文件的方法;
2、开始养成给代码做注释的习惯;
3、通过视频学习了写可复用代码的思路。
回复 支持 反对

使用道具 举报

0

主题

4

帖子

164

积分

注册会员

Rank: 2

积分
164
发表于 2017-2-12 16:39:36 | 显示全部楼层
[Python] 纯文本查看 复制代码
import glob

my_diction ={} #写字典表

def Readfile(file):
    file_obj = open(file)
    try:
        while True:
            ##读取文件并分割
            line = file_obj.readline().strip("\n")#按照行来读取文件并进行去换行符
            if not line:
                break                             #如果文件行读取完,跳出读取文件
            array = line.split("\t")              #分割键值
            ##对相同的键进行扩展,存入my_diciton内,返回至全局;对未出现的基因添加新的键值对
            if array[0] in my_diction:
                my_diction[array[0]].append(array[1])
            else:
                my_diction[array[0]] = [array[1]]
    finally:
        file_obj.close()
    return my_diction # 返回值


def main():
    list_dirs = glob.glob("./*.txt") #得到当前目录下的txt文件并存入变量list_dirs
    for i in list_dirs:
        Readfile(i)
    with open('./sample.csv','w') as fp_out:#准备当前目录下sample.csv文件的句柄
        for gene_name in my_diction:
            print("%s\t%s" % (gene_name,"\t".join(my_diction[gene_name])))
            fp_out.write("{},{}\n".format(gene_name, ",".join(my_diction[gene_name])))#写入csv

if __name__ == '__main__':
    main() #调用函数

小结:1 字典表在全局变量,便于利用循环对值进行扩增
         2 读取文件的选择方式
         3 只有反复敲代码!!才有顿悟!!
回复 支持 反对

使用道具 举报

2

主题

21

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-2-12 20:37:28 | 显示全部楼层
本帖最后由 AdaWon 于 2017-2-12 21:18 编辑

147-python
1. 首先使用htseq-count分开对每个样本计算所有的基因表达量,生成一个个独立的文件;其次,第一列是基因,第二列是该基因的counts值;
如:
$ ls
0dpa_R1.count   1dpa_R2.count  2dpa_R3.count  3dpa_R1.count  4dpa_R2.count   
0dpa_R2.count   1dpa_R3.count  3dpa_R2.count  4dpa_R3.count  5dpa_R3.count
0dpa_R3.count   2dpa_R1.count  3dpa_R3.count  5dpa_R1.count   
1dpa_R1.count  2dpa_R2.count   4dpa_R1.count  5dpa_R2.count


$ cat -A 0dpa_R1.count |more
gene:Solyc00g005000.2^I0$
gene:Solyc00g005020.1^I0$
gene:Solyc00g005040.2^I0$
gene:Solyc00g005050.2^I0$
gene:Solyc00g005060.1^I0$
gene:Solyc00g005070.1^I0$
gene:Solyc00g005080.1^I0$
gene:Solyc00g005090.1^I0$
gene:Solyc00g005100.1^I0$
gene:Solyc00g005110.1^I0$
gene:Solyc00g005120.1^I0$
gene:Solyc00g005130.1^I0$
gene:Solyc00g005140.1^I0$
gene:Solyc00g005150.1^I0$



$  less -S 0dpa_R1.count
gene:Solyc01g006780.2   9
gene:Solyc01g006790.2   2
gene:Solyc01g006800.2   6
gene:Solyc01g006810.2   0
gene:Solyc01g006820.2   8
gene:Solyc01g006830.2   4
gene:Solyc01g006840.2   10
gene:Solyc01g006850.1   0
gene:Solyc01g006860.2   0
gene:Solyc01g006870.2   0
gene:Solyc01g006880.2   0
gene:Solyc01g006890.2   3
gene:Solyc01g006900.2   6
gene:Solyc01g006910.2   0
gene:Solyc01g006920.2   0
gene:Solyc01g006930.2   0
gene:Solyc01g006940.2   20
gene:Solyc01g006950.2   1
gene:Solyc01g006960.2   4
gene:Solyc01g006970.2   2
gene:Solyc01g006980.2   240#因私人保密性,更改了部分counts值

2.将多个文件合并成一个大的行列式,基因名作为row name, 每个文件的counts 值对应的文件作为col name;
代码如下:

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

import glob #

MyDiction = {} 
#设定空字典,字典是一种可变容器模型,且可存储任意类型对象;键必须是唯一的,但值则不必;值可以取任何数据类型,但键必须是不可变的,如字符串,数字或元组。

def Readfile(file): #定义 Readfile 函数,读取文件,将其中信息存入MyDiction
        file_obj = open(file)
        try:
                while True:
                        line = file_obj.readline().strip('\n')  #读一行操作一行,即对每一行去除换行符
                        if not line:  #如果此行为空,则跳出循环 or len(line) == 0:
                                break
                        array = line.split('\t')  #对每一行进行分割操作,以空格键作为分隔符,返回分割后的字符串形成一个数组array
#Python的数组分三种类型:list_arr = [元素],tuple_arr = (元素),Dictionary_arr = {元素k:v}
                        
                        if array[0] in MyDiction:  #array[0] 为基因名,即为键值key
                                MyDiction[array[0]].append(array[1])   #array[1]为counts值,即为value
                        else:
                                MyDiction[array[0]] = [array[1]] 
        finally:
                file_obj.close()
        return MyDiction

def main(): #构建主函数体
        list_dirs = glob.glob('./*.count') #all .count files in pwd
#glob包最常用的方法只有一个, glob.glob();功能与Linux中的ls相似,用于查询目录下的文件,即列出所有符合该表达式的文件(与正则表达式类似),将所有文件名放在一个表中返回。
        for i in list_dirs:
                Readfile(i)
        for key in MyDiction.keys(): #遍历键值
                print "%s\t%s"%(key,'\t'.join(MyDiction[key])) #格式化输出 ""%(); join 连接


if __name__ == '__main__':
        main()


3. 运行
$ python fourth_task.py  > all_htseq_cout_result.txt
$ less -S all_htseq_cout_result.txt
......
gene:Solyc03g093730.1   0       0       0       0       0      880    0       0       0       0       0       0 ...
gene:Solyc12g009710.1   0       90     0       0       0       0       0       0       0       0       0       0 ...
gene:Solyc03g046170.1   0       0       0       0       0       0       47     0       0       0       10     0 ...
gene:Solyc06g009710.2   4       0       0       0       0       0       0       0       0       0       1       0 ...
gene:Solyc02g069560.2   0       0       0       22     12     0       0       0       0       0       0       6 ...

......

回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-2-13 20:12:53 | 显示全部楼层
本帖最后由 Aiyawq 于 2017-2-14 10:16 编辑

207-perl-wangQ;
哎呀,不得不说二维哈希真的是好用,上一节课用了,这一回居然还发现这么好一个功能;代码如下,给各位分享一下吧~

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my $indir = $ARGV[0];
opendir DIR,$indir;
my (%hash1,%hash2,$c);
open A,">result.txt";
foreach my $file (sort {$a cmp $b} grep(/^\S+\.txt$/,readdir(DIR))){
        open IN,"$indir/$file";
        while (<IN>){
                chomp;
                my @c=split(/\t/,$_);
                $hash1{$c[0]}{$file}=$c[1];
                $hash2{$file}=1;
        }
}
my @d=sort keys %hash2;
my $d=join "\t",@d;
print A "EnsEMBL_Gene_ID\t$d\n";
close $indir;
closedir DIR;
foreach my $e(sort keys %hash1){
        print A $e;
        foreach my $f(sort keys %{$hash1{$e}}){
                print A "\t$hash1{$e}{$f}";
        }
        print A "\n";
}
__END__



本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

5

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2017-2-14 10:32:39 | 显示全部楼层
109-R+perl+python-死小孩
[Python] 纯文本查看 复制代码
#coding=utf-8
import os
import os.path
from _collections import OrderedDict
import re

print("使用说明:\n\t输入需要合并文本的文件夹的路径,再输入每个文本A1所不需要的序列,在对应的路径产生目标文本")

sample=OrderedDict()
genes=OrderedDict()

path=input('请输入路径')
title=input('请输入标签')

path=re.sub(r'\\','\\\\\\\\',path)
g=open(path+'\\express.txt','w')
files=os.listdir(path)#得到path路径下的所有文件名称
for file in files:
    with open(path+'/'+file,'r') as gene:
        for line in gene:
            line=line.rstrip()
            ge=line.split('\t')[0]
            val=line.split('\t')[1]
            if ge==title:#读取第一行,并以sample+样品名建立字典,字典key为基因名,value为表达量
                sam=val
                sample[sam]=1
                sample[sam]=OrderedDict()
            else:#数据写入字典
                sample[sam][ge]=val
                if ge not in genes.keys():  
                    genes[ge]=1
g.write("\t")  #输出第一行标题      
for name in sample:
    g.write(name)
    g.write("\t")
g.write('\n')
for k in genes:#输出每行数据
    li=k+'\t'
    for l in sample:
        li=li+str(sample[l][k])+'\t'
    li=li+'\n'
#    print(li)   
    g.write(li)        
g.close()
回复 支持 反对

使用道具 举报

0

主题

6

帖子

337

积分

中级会员

Rank: 3Rank: 3

积分
337
发表于 2017-2-14 13:14:51 | 显示全部楼层
猪兔子2号 发表于 2017-1-4 22:44
htseq-count出来一样的行的,我会用下面这个命令来做
[mw_shl_code=shell,true]paste file1 file2 file3 |  ...

楼主,最后一个 print $i 也可以写成 print ""。这一行处理之后,输出换行符,再处理下一行。
回复 支持 反对

使用道具 举报

0

主题

6

帖子

337

积分

中级会员

Rank: 3Rank: 3

积分
337
发表于 2017-2-14 13:41:46 | 显示全部楼层
035-python-shell
我之前在做RNA-seq的时候,用的是 DESeq2这个包。在用这个包的时候,里面有个DESeqDataSetFromHTSeqCount 这个命令,可以将这些文件导入。看到大家的
答案,获益很多。python在努力的学习当中,目前还写不出来。
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.