搜索
查看: 7157|回复: 79

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

  [复制链接]

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
发表于 2016-12-29 09:19:19 | 显示全部楼层 |阅读模式
相信用过htseq-count的朋友都知道,它是分开对每个样本计算所有的基因表达量,所以会生成一个个独立的文件,我用perl脚本模仿它的结果如下:
$ head a.txt
gene_1  178
gene_2  692
gene_3  486
gene_4  666
gene_5  395
gene_6  48
gene_7  926
gene_8  733
gene_9  660
gene_10 578

第一列是基因,第二列是该基因的counts值,共有a~z这26个样本的counts文件,需要合并成一个大的行列式,这样才能导入到R里面做差异分析,如果手工用excel表格做,当然是可以的,但是太麻烦,如果有500个样本,正常人都不会去手工做了,需要编程。

生成测试文件的代码如下:
[AppleScript] 纯文本查看 复制代码
#首先新建文件tmp.sh 输入这个代码:
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
## 然后用perl脚本调用这个tmp.sh文件:
perl -e 'system(" bash tmp.sh >$_.txt") foreach a..z'
##这样就生成了a~z这26个样本的counts文件


用shell或者perl或者python,设置R语言都可以做,但是各有优缺点,而且如果每个样本的基因顺序并不一致,这时候你应该怎么做呢?

实际需求如下:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213 里面有56个文件(ftp://ftp.ncbi.nlm.nih.gov/geo/s ... pl/GSE48213_RAW.tar),需要合并成一个表达矩阵,来根据cell-line的不同,分组做差异分析。
paper是:https://www.ncbi.nlm.nih.gov/pubmed/24176112

输出的表达矩阵,如下所示:


请务必做出此题,很重要!


先给一下shell结合R语言的做法:
[Python] 纯文本查看 复制代码
## 首先在GSE48213_RAW目录里面生成tmp.txt文件
awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt
## 然后把tmp.txt导入R语言里面用reshape2处理即可!
setwd('tmp/GSE48213_RAW/')
a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)



本帖子中包含更多资源

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

x



上一篇:看了一篇HATs和HDACs的文章
下一篇:结构变异SV是什么?
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

95

主题

119

帖子

614

积分

版主

Rank: 7Rank: 7Rank: 7

积分
614
发表于 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

使用道具 举报

0

主题

3

帖子

28

积分

新手上路

Rank: 1

积分
28
发表于 2017-1-6 19:06:01 | 显示全部楼层
I think pandas can make our job easier!

[Python] 纯文本查看 复制代码
import pandas
import os
name_list=os.listdir("GSE48213_RAW")
fram_list=[pandas.read_table("GSE48213_RAW/%s"%name) for name in name_list]
fram=fram_list[0]
for i in range(1,len(fram_list)):
    fram=pandas.merge(fram,fram_list[i])
fram.to_csv("result.csv",index=False)


by 075-python-哈哈哈
回复 支持 3 反对 0

使用道具 举报

1

主题

2

帖子

34

积分

新手上路

Rank: 1

积分
34
发表于 2016-12-29 11:12:00 | 显示全部楼层
本帖最后由 yydBIG 于 2016-12-29 11:13 编辑

贴一个perl实现的版本,不用排序,没出现的用0代替。使用方法:
perl zuhe.pl file1 file2 ...
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
#2010-12-24
#yangyd

my %total;
my $header = "Gene";

foreach my $i(0..@ARGV-1){
  $header.="\t$ARGV[$i]";
  open I,"$ARGV[$i]";
  while(<I>){
    chomp;
    my @tmp=split;
    $total{$tmp[0]}[$i]=$tmp[1];
  }
}
print $header,"\n";
foreach my $key(keys %total){
  foreach my $i(0..@ARGV-1){
    if(!$total{$key}[$i]){
      $total{$key}[$i]=0;
    }
  }
  print $key,"\t",join("\t",@{$total{$key}}),"\n";
}


回复 支持 2 反对 0

使用道具 举报

5

主题

36

帖子

435

积分

中级会员

Rank: 3Rank: 3

积分
435
发表于 2017-2-28 17:08:46 | 显示全部楼层
本帖最后由 W.Peng 于 2017-2-28 17:12 编辑

合并文件可以用linux的join命令,将两个文件里指定栏位同样的行连接起来,但这两个文件必须是按照同样的规则排序的,空值可以用指定的字符(如 NA )代替。
[Bash shell] 纯文本查看 复制代码
#!/usr/bin/bash

for i in /data_dir/*.txt
do
    sort ${i} > ${i}.new
    awk '{print$1}' ${i}.new > tmp.txt
done

for j in /data_dir/*.txt.new
do
    join -e 'NA' tmp.txt ${j} > join.txt
    mv join.txt tmp.txt
done
mv tmp.txt result.txt ; rm -rf *.txt.new
回复 支持 1 反对 0

使用道具 举报

1

主题

16

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-8 20:13:57 | 显示全部楼层
本帖最后由 disheng 于 2017-1-8 20:37 编辑

043-perl-R-shell-涤生

学习了二楼yydBIG的二维哈希,真的太棒了。参照他的思想,用哈希也写了个
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

my %hash;
my $header = "Gene";
foreach my $i ( 0..@ARGV-1){
        $header .= "\t$ARGV[$i]";
        open I,"$ARGV[$i]";
        while(<I>){
                chomp;
                my($gene_num,$num) = split;
                $hash{$gene_num} .= "$num\t";
        }
}
print $header,"\n";
foreach $i (keys %hash){
        print "$i\t$hash{$i}\n";
}


本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

0

主题

10

帖子

101

积分

注册会员

Rank: 2

积分
101
发表于 2016-12-29 09:32:23 | 显示全部楼层
大多数情况下基因顺序是一致的
回复 支持 反对

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 2016-12-29 09:41:44 | 显示全部楼层
liyan 发表于 2016-12-29 09:32
大多数情况下基因顺序是一致的

这种情况下是一致的,但是,如果呢,比如其它软件,其它情况下,如果你需要合并vcf文件,坐标就不一致了,怎么办,能hold住吗?
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

9

帖子

90

积分

注册会员

Rank: 2

积分
90
发表于 2016-12-29 10:26:34 | 显示全部楼层
本帖最后由 l0o0 于 2016-12-29 10:49 编辑

这里我觉得用R语言里面的merge来做感觉会比perl或python建字典方便许多。
这里不方便贴代码,可以看下github上的链接https://github.com/l0o0/bio-analysis-kit/blob/master/merge.R,默认文件第一行是表头。
[AppleScript] 纯文本查看 复制代码
#!/bin/env Rscript

args = commandArgs(T)

if (length(args) != 4){
    print('USAGE: Rscript merge.R pattern merge_colname input_dir output_file')
    quit()
}


files = list.files(args[3],pattern = paste(args[1],'$',sep=''))

a <- files[1]
atable <- read.table(a, header=T, check.names=F, sep='\t')
samplea = strsplit(basename(a), '.', fixed=T)[[1]][1]
colnames(atable) = c(args[2], samplea)

for (i in files[-1]){
	print (i)
	sampleb = strsplit(basename(i), '.', fixed=T)[[1]][1]
    btable <- read.table(i, header=T, check.names=F, sep='\t')
    colnames(btable) = c(args[2], sampleb)
	print(colnames(btable))
	#names(btable) <- c('gene_id', fpkm)
	atable <- merge(atable, btable, by = args[2], all =TRUE)
}

row.names(atable) = atable[,1]
atable = atable[,-1]
print(colnames(atable))
write.table(atable, sep='\t', file=args[4], col.names=NA,quote=FALSE, na='NA')

Rscript merge.R file-pattern  geneid $PWD test.txt

回复 支持 反对

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 2016-12-29 10:31:33 | 显示全部楼层
l0o0 发表于 2016-12-29 10:26
这里我觉得用R语言里面的merge来做感觉会比perl或python建字典方便许多。[mw_shl_code=applescript,true]ar ...

我给你示范了一下如何排版,你继续修改一下你的代码,稍微正常一点,编辑器里面有插入代码选项
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 2016-12-29 15:46:05 | 显示全部楼层
yydBIG 发表于 2016-12-29 11:12
贴一个perl实现的版本,不用排序,没出现的用0代替。使用方法:
perl zuhe.pl file1 file2 ...
[mw_shl_cod ...

你的代码很简练,赞,虽然我没有测试过,但是跟我的思路是一样的
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

3

主题

10

帖子

94

积分

注册会员

Rank: 2

积分
94
发表于 2017-1-4 22:44:56 | 显示全部楼层
htseq-count出来一样的行的,我会用下面这个命令来做
[Shell] 纯文本查看 复制代码
paste file1 file2 file3 | awk '{printf $1"\t";for(i=2;i<=NF;i=i+2)printf $i"\t";print $i}'
回复 支持 反对

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

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

你这个也很赞
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

41

帖子

266

积分

中级会员

Rank: 3Rank: 3

积分
266
发表于 2017-1-6 10:54:28 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-6 11:12 编辑

059pyth+R 冷洋
意思是要把所有的文本文件的内容合成一个文本文件吗?
[Python] 纯文本查看 复制代码
#coding=utf-8
import os
import os.path
g=open('D:/one.txt','w')
path='D:/56-breast-cancer-celline-RNA-seq'
files=os.listdir(path)#得到path路径下的所有文件名称
for file in files:
    f=open(path+'/'+file)#逐个打开文件
    for i in f:
        g.write(i+'\n')#将每个文件的内容写入one.txt
g.close()

回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-10-17 09:56 , Processed in 0.040957 second(s), 43 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.